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Constraint Treatment Techniques and Parallel Algorithms 
for Multibody Dynamic Analysis 
Dissertation directed by Professor K. C. Park 

Computational procedures for kinematic and dynamic analysis of 
three-dimensional multibody dynamic (MBD) systems are developed from 
the differential-algebraic equations (DAEs) viewpoint. First, to minimize 
constraint violations during the time integration process and to obviate de- 
graded constraint force solution involving ill-conditioned matrices, two ro- 
bust and efficient constraint treatment techniques, t. e., penalty constraint 
stabilization technique and natural partitioning scheme, are developed. The 
computational issues for enhancing accuracy, stability, and programming 
modularity of the techniques are also addressed for MBD analysis. 

Second, to treat the governing equations of motion, a two-stage 
staggered explicit-implicit numerical algorithm, that takes advantage of a 
partitioned solution procedure and a robust and parallelizable integration 
algorithm, is developed. Mainly, this algorithm uses a two-stage staggered 
central difference algorithm to integrate the translational coordinates and 
the angular velocities. The angular orientations of bodies in MBD systems 
are then obtained by using an implicit algorithm via the kinematic rela- 
tionship between Euler parameters and angular velocities. It is shown that 
the combination of the present solution procedures yields a computationally 
more accurate solution. 

Third, to speed up the computational procedures, parallel imple- 
mentation of the present constraint treatment techniques, two-stage stag- 
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gered explicit-implicit numerical algorithm has been efficiently carried out. 
To this end, the DAEs and the constraint treatment techniques have been 
transformed into arrowhead matrices to which Schur complement form has 
been derived. By fully exploiting the sparse matrix structural analysis tech- 
niques, a parallel preconditioned conjugate gradient numerical algorithm is 
used to solve the systems equations written in Schur complement form. 

To evaluate the computational procedures developed in the present 
work, a software testbed has been designed and implemented in both sequen- 
tial and parallel computers. This testbed has been used to demonstrate the 
robustness and efficiency of the constraint treatment techniques, the accu- 
racy of the two-stage staggered explicit-implicit numerical algorithm, and 
the speed up of the Schur-complement-based parallel preconditioned conju- 
gate gradient algorithm on a parallel computer. 
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CHAPTER I 


INTRODUCTION 


1.1 OVERVIEW 

The kinematic and dynamic analysis of three dimensional mnltibody 
dynamic (MBD) systems has attracted the attention of many researchers 
over the past two decades. This is due to the fact that many mechanical 
systems of interest in industry can be effectively modeled by using systems 
of linked bodies. Moreover the rapid development of computer hardware 
and software has also played an important role in making the computer 
simulation of MBD systems more realistic if the number of bodies in the 
systems remain small. These research activities have primarily concentrated 
on improving either the design and verification of the control system, or the 
system design and dynamic analysis of multidisciplinary engineering prob- 
lems. As a result, several stand-alone general-purpose computer programs 
[1-11] which are based on different approaches have been developed. 1 hose 
computer programs possess the capability to automatically generate and nu- 
merically integrate the equations of motion of multibody problems such as 
robot arm maneuvers, spacecraft dynamics, and ground vehicle dynamics. 

However, when these systems become complex, computational effi- 
ciency becomes a dominant issue during the preliminary design stage that 
may require many analysis iterations. This has motivated several research 
groups to make effective use of parallel computational technology 12-11; in 
order to speed up the dynamics analysis of MBD systems, thus ultimately 
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achieving real-time simulation for large-scale problems. The key issues in 
exploiting the parallelism inherent in MBD systems include a versatile data 
structure for describing system topology, an automatic procedure to gen- 
erate the system equations of motion, a streamlined treatment of system 
constraints, a robust time integration algorithm, and an easy interpretation 
of the simulation results. 

1.1.1 Systematic Study of MBD Formulations 

In general, the equations of motion for MBD systems can be derived 
and expressed in various forms depending upon the type of coordinates 
chosen to describe the configuration of the bodies. An important kinematic 
characteristic of these coordinates is how they treat the joints that, are used 
to describe the kinematic relationships of the bodies in the systems. 1 hus if 
an arbitrary set of coordinates is chosen, the final system dynamic equations 
can be interpreted as results of two basic approaches: the augmentation 
approach and the elimination approach as shown in Fig. 1.1. The first 
approach gives a set of differential-algebraic equations whereas the second 
approach gives a set of second-order differential equations. 

In order to understand the advantages and disadvantages of using 
different coordinates to derive the equations of motion, four choices ot co- 
ordinates will be discussed. The first choice is to use a set of independent 
coordinates , which determine the position of bodies with the least possi- 
ble number of state variables. A minimal set of second-order differential 
equations, which is given in terms of system independent variables, is ob- 
tained in which the constraint conditions are absent. However, the rapidly 
growing complexity in the derivation as the number of variables increases. 
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Fig. 1.1 MBD Equations of Motion and Their Corresponding 
Computational Procedures 
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and the high degree of nonlinearity of the equations of motion make these 
coordinates difficult to implement in a general-purpose computer program. 

The second choice is that of relative coordinates [1-4,15,16], which 
define the orientation of each moving body with respect to either a non- 
moving body or another adjacent moving body. For an open tree structure, 
the number of relative coordinates is equal to the number of independent 
coordinates. For a closed-loop system, constraint equations are imposed 
via Lagrange multipliers, in which case the number of relative coordinates 
exceeds that of independent coordinates. Relative coordinates have the 
disadvantage in that they do not directly determine the position of each 
body in the system; thus postprocessing of the simulation results is needed. 
Furthermore, when the system consists of several closed loops, extensive 
preprocessing is needed to identify an appropriate set of independent vari- 
ables. 

The third choice is that of natural coordinates [17.18], which define 
a body using two or more moving coordinates rigidly attached to it. 1 hose 
moving coordinates are located preferably at the joints of the mechanism, 
and can be shared by adjacent bodies. The main advantage of natural co- 
ordinates is that they lead to a simple computer implementation and easy 
formulation in conjunction with quadratic or linear constraint equations. 
However, the presence of a fully populated mass matrix renders these coor- 
dinates less attractive in parallel computation. Another drawback ol these 
coordinates is that during the process of numerical integration a position 
can be reached which causes the matrix that is used to identify the variable 
dependencies to become singular. Should this happen, a new linear combi- 
nation matrix need to be constructed in order to continue the integration 
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process. 

The last choice is that of Cartesian coordinates [5-8], which define 
the position of each particle in each individual body in the system with re- 
spect to an inertial reference frame. The angular orientation of each body is 
defined by the body-fixed reference frame via Euler parameters or Euler an- 
gles. The main advantages of this choice is that the equations of motion are 
easy to derive, which facilitates the development of general-purpose com- 
puter programs. Since these coordinates yield a maximal set of equations, 
redundant coordinates and Lagrange multipliers have to be solved as part of 
the simulation process, which may lead to computational inefficiency unless 
special attention is paid to computational issues. 

If independent coordinates are used, the equations of motion are 
generated in terms of system degrees of freedom expressed in differential 
equation form. Obviously this approach leads to a minimal set of equations 
of motion but suffers from the appearance of dense solution matrices and 
highly nonlinear kinematic descriptions. 

When d’Alembert’s principle of virtual work together with Lagrange 
multipliers are applied to the systems based on relative coordinates, natu- 
ral coordinates or Cartesian coordinates, the resulting equations for MBD 
systems are given by a set of second-order differential equations augmented 
with algebraic constraint equations. These combined system of equations 
belong to the class of differential-algebraic equations (DAEs). In contrast 
to the independent coordinates approach, DAEs make use of a larger num- 
ber of equations yet preserve the sparsity of the solution matrix as well as 
the simplicity of the kinematic relationships. Furthermore, the approach is 
amenable to implementation in modular, general-purpose MBD programs. 
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1.1.2 Evaluations of Existing Computational Procedures 

A closed-form solution of the MBD equations is in general impossi- 
ble except for highly simplified problems. Thus time integration algorithms 
must be used to obtain the numerical solution of the system governing equa- 
tions. For the second-order differential equations produced by the elimina- 
tion approach, both the modified explicit central difference formula and 
as well as stiffly-stable formulas in conjunction with the Newton-Raphson 
algorithm may yield reasonably stable and accurate solutions. As for dif- 
ferential algebraic equations, Gear [19,20] has investigated a special class 
of numerical algorithms for the solution of some restricted DAE problems. 
Orlandea et al. [6] have applied this solution technique together with a 
sparse matrix formulation but encountered numerical problems because the 
discrete system equations to be solved often become numerically stiff and 
ill-conditioned. 

An alternative approach, advocated by Gear and Petzold j21,22i, 
relies on augmenting the second-order governing equilibrium equations with 
twice time-differentiated constraint equations so that numerical ordinary dif- 
ferential equations solvers can be applied. However, numerical integration 
algorithms provide only an approximate solution. As a result, numerical 
errors will propagate and accumulate so that eventually the constraint con- 
ditions are no longer satisfied within the desired accuracy. One approach to 
stabilize the constraint violations was proposed by Baumgarte [23,24] who 
modified the original constraint equations to form a set of relaxation dif- 
ferential constraint equations. Park and Chiou [25,26] have shown that for 
some MBD problems Baumgarte’s constraint stabilization technique suflers 
from ill-conditioning in the solution for Lagrange multipliers. 1‘ urthei more. 
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for complicated MBD systems, the process of determining optimal relax- 
ation parameters that are used to tailor the constraint violations to each 
specific problems may encounter computational difficulties. 

An ultimately different approach to avoid constraint violations con- 
sists of eliminating the Lagrange multipliers from DAEs so that a set of 
second-order differential equations is obtained. This can be done by identi- 
fying system dependent and independent variables from the given constraint 
Jacobian matrix so that the null space of the constraint Jacobian matrix 
can be formed and consequently used to eliminate the Lagrange multipliers. 
In order to find a set of numerically superior independent variables, sev- 
eral numerical algorithms have been employed to decompose the constraint 
Jacobian matrix. These algorithms include: the generalized coordinate par- 
titioning scheme [27], the singular value decomposition [28,29k the natural 
coordinates partitioning scheme [17,18], and the null space scheme [30-32]. 
Since these computational schemes for determining the set of independent 
coordinates can become computationally expensive, the chosen set of inde- 
pendent coordinates is maintained during the numerical simulation until the 
specified accuracy criteria are violated. When this occurs, it is necessary to 
choose a new set of independent coordinates by repeating the idem ilicat ion 
process. 

Recently, methods based on 0{n ) algorithms, where n is the number 
of generalized coordinates, and several variations have been proposed [33- 
39]. These algorithms are primarily applicable to MBD systems consisting of 
tree topologies in which their equations of motion may be recursively solved 
in O(n) operations. If the system topology embodies multiple closed loops, 
significant modifications are required in order to obtain numerical solutions. 
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Moreover, the presence of closed loops may cause O(n) algorithms to loose 
the simplicity of open tree topology in parallel computations. 

1.2 OBJECTIVES 

Since the aforementioned solution procedures suffer various draw- 
backs in the computer implementation, these have motivated us to look for 
alternative solution procedures that overcome those difficulties. Alternative 
solution procedures that involve either constraint stabilization or constraint 
elimination overcome the following disadvantages: unacceptable constraint 
violation during the process of time integration; degraded constraint force 
solution involving ill-conditioned matrix; large computational expense in 
computing the null space of the constraint Jacobian matrix; inefficiency in 
using the implicit iterative algorithms; and difficulties in extending existing 
algorithms to parallel computations. 

The objectives of this dissertation are to develop: first, robust and 
efficient treatments of constraints; second, explicit-implicit integration algo- 
rithms to solve DAEs efficiently and accurately; and third, a parallel imple- 
mentation procedure for a general multibody dynamics simulation capacity. 

With these objectives in mind, we will first review the spatial kine- 
matics of linked bodies by employing two sets of coordinates that describe 
the configuration and velocity of bodies in the system. Inertial coordinates 
are adopted to locate the center of mass of each of the bodies. Body-fixed 
coordinates are rigidly attached to the center of mass of each body so that 
angular orientations of the bodies in the system can be obtained as soon as 
angular representations are determined. Note that the purpose ot chosing 
inertial coordinates for the translational motions and body-fixed coordi- 
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nates for the angular motions is to decouple the inertia matrix to obtain the 
translational and rotational equations separately. After the kinematics of 
a body in space has been determined, a formulation based on d’Alembert’s 
principle of virtual work is adopted to derive the system governing equa- 
tions by treating the bodies in the system as originally independent of each 
other. Nonlinear kinematic constraint conditions, which are imposed to de- 
scribe the interconnectivity between the various bodies, are appended to 
the formulation using Lagrange multipliers. The final system equations of 
motion, which are characterized a s DAEs, not only enhance programming 
modularity but can also be generated automatically. As mentioned previ- 
ously, the use of existing numerical time integration solution procedures may 
encounter computational difficulties. In this regard, two newly developed 
schemes based on constraint stabilization (penalty constraint stabilization 
scheme) and constraint elimination (natural partitioning scheme), are intro- 
duced to correct for the constraint violations accurately and efficiently. 

1.2.1 Robust and Efficient Treatment of Constraints 

The penalty constraint stabilization scheme is based on the obser- 
vation that time-differentiated equations of penalty form retain a parabolic 
characteristic in time. Thus, as time progresses constraint violations will de- 
cay according to intrinsic time constants. This penalty time-differentiated 
form, which is given by the time rate of the constraint forces, enables us to 
overcome the difficulties that have been encountered in Baumgarte’s tech- 
nique. Furthermore, this scheme offers the attractive feature that the sys- 
tem equations can be processed in two modules: the generalized coordinates 
module and the generalized constraint forces module. This separation fits 
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nicely in the framework of the partitioned procedure [40-42] adopted for the 
time integration. 

The natural partitioning scheme, which is quite different from the 
penalty constraint stabilization scheme, uses the existing physical coordi- 
nates by explicitly identifying their dependent and independent coordinates 
without relying heavily on the numerical algorithms that have been men- 
tioned previously. The resulting null space of the constraint Jacobian matrix 
can be generated in parallel if the system topology consists of several open 
chains. Applying the null space of the constraint Jacobian matrix to the gov- 
erning DAEs leads to the elimination of the Lagrange multipliers and yields 
a set of second-order differential equations that are expressed in terms of 
system independent variables. 

1.2.2 Explicit-Implicit Integration of MBD Equations 

In the present formulation, both angular velocity-dependent cen- 
tripetal accelerations and the angular accelerations appear in the equations 
of motion that represent the rotational motions of linked bodies. Direct time 
integration of angular velocities, except for some simple kinematic configu- 
rations, does not directly yield the angular orientation of a body in space. 
Hence, a partitioned solution procedure [40-42], which has the capability 
to separately solve coupled systems of equations while treating interaction 
terms as external forces, is used to separately integrate translational and 
rotational equations. To obtain a robust and parallelizable integration algo- 
rithm, the explicit central difference time integration formula is thought to 
be the best candidate to treat the partitioned translational and rotational 
quantities. If the central difference formula is adopted, the approximation 
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of the angular velocity for the evaluation of angular acceleration leads to 
numerical instability for the governing equations of motion. This has moti- 
vated us to exploit a two-stage explicit procedure which stabilizes the central 
difference algorithm and delivers an accurate solution. The Euler param- 
eters are used in the present solution procedure to represent the angular 
orientations of bodies. An implicit mid-point time integration algorithm 
is employed to integrate the Euler parameters by exploiting the kinematic 
relationship with the associated angular velocities. The specific implicit 
algorithm presented has been chosen because it is unconditionally numer- 
ically stable while it can be analytically inverted during actual computer 
implementation because of its special four by four matrix form. 

Combining these solution algorithms a two-stage staggered explicit- 
implicit solution procedure [43] has been developed. This solution proce- 
dure, which invokes either the penalty constraint stabilization scheme or 
the natural partitioning scheme to stabilize the constraint violations, has 
been implemented in a computer program to validate and demonstrate its 
robustness and accuracy. 

The present solution procedure based on the penalty constraint sta- 
bilization scheme consists of two modularized solvers: the generalized coor- 
dinate solver and the constraint force solver. The solution procedure based 
on the natural partitioning scheme includes the generalized coordinate solver 
and the independent coordinate solver. The generalized coordinate solver 
combines an improved version of the explicit central difference algorithm 
for integration of the translational coordinates and angular velocity with an 
implicit algorithm to update the Euler parameter representation of angular 
orientations by exploiting the uncoupled inertia expression. The procedure 
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has successfully been interfaced with the penalty staggered stabilization 
technique which solves the constraint forces as independent variables by 
implicitly integrating a stabilized companion differential equation for the 
constraint forces in time. The combination of the two algorithms can be 
invoked in a sequential manner on the rigid and flexible components of the 
multibody system resulting in an attractive, modular solution procedure. 
As for the independent variable solver, a procedure based on body-by-body 
constraint Jacobian matrices is developed to explicitly form the null space of 
the constraint Jacobian matrix and consequently obtain system independent 
variables. 

1.2.3 Parallel Implementation Procedures for MBD Analyzer 

Since an MBD system may consist of hundreds or even thousands of 
bodies, the numerical solution may consume a prohibitive amount of CPU 
time. To reduce the CPU time dramatically, it is advantageous to develop 
efficient parallel algorithms by incorporating existing parallel computers. In 
general, issues involving parallel computations of MBD systems include gen- 
eration of the system equations of motion, incorporation or elimination of 
constraint forces, integration of generalized coordinates, and interpretation 
of the simulation results. A method for exploiting the parallelism of the 
present constraint stabilization, constraint elimination and two-stage stag- 
gered explicit-implicit solution procedure has been developed. The main 
thrust of this method uses a Schur-complement-based parallel precondi- 
tioned conjugate gradient numerical algorithm to decide either the gen- 
eralized acceleration vector and constraint forces of the penalty constraint 
stabilization scheme or the generalized and independent acceleration vec- 
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tors of the natural partitioning scheme. The present algorithm has been 
implemented and tested on existing parallel computers and has provided 
encouraging results in practical MBD problems. 

1.3 DISSERTATION OUTLINE 

The dissertation is organized as follows. Chapter 2 reviews the spa- 
tial kinematics of rigid bodies used in the present MBD formulation. Chap- 
ter 3 employs d’Alembert’s principle of virtual work to derive the governing 
equations of motion that consist of a set of algebraic constraint equations 
coupled with the second-order differential equations of motion. Chapter 4 
derives the mechanical properties of the joints that connect bodies in the 
MBD system. These joints are introduced in the dynamic formulation using 
a set of algebraic constraint equations that are adjoined to the equations of 
motion in constraint Jacobian matrix form. Chapter 5 deals with the DAEs 
by reviewing several existing numerical solution procedures. Two newly de- 
veloped schemes based on the penalty constraint stabilization scheme and 
the natural partitioning scheme are employed to overcome computational 
difficulties associated with those solution procedures. A two-stage staggered 
explicit-implicit algorithm for updating the translational coordinates, angu- 
lar orientations, and constraint forces is developed based on the proposed 
technique. Chapter 6 analyzes and exploits the parallelism inherent in the 
solution procedures. Chapter 7 gives numerical example problems in order 
to demonstrate the robustness and efficiency of utilizing the present, solu- 
tion procedures. Chapter 8 summarizes the accomplishments of the present 
investigation and discusses directions for further research in the field of 
multibody systems. 





CHAPTER II 


SPATIAL KINEMATICS 


2.1 Introduction 

Kinematics, which is the study of the motion of particles and bodies 
without the forces associated with these motions, has been used to analyze 
the position, velocity, and acceleration of bodies and determine the design 
geometry of the bodies in the mechanical systems. In this chapter, we 
begin with the derivation of different finite rotational representations and 
subsequently obtain the position, velocity, and acceleration of a particle in 
space. Finally, we will consider the particle as if it has been attached to 
a rigid body and thus ultimately complete the derivation of kinematics for 
the rigid body. 

2.2 Reference Frames 

In mechanics, a most fundamental technique is using vectorial quan- 
tities to locate the position of a particle in a given reference frame. When 
the position vector from the origin of that reference frame to the particle 
has been defined, we can resolve this position vector by one or more systems 
of coordinates for a particular use. In many dynamics problems, relations 
between the component of the vector in various reference frames prove ex- 
tremely useful. To derive such relations, let us consider a position vector r., 
(Fig. 2.1) expressed in terms of three components parallel to the three axes 
of a Cartesian frame: 


r P = + r e 2 e 2 + r 3 e 3 


( 2 . 2 . 1 ) 



psaor*; 


t-LANK NOT FILMED 


15 

where (ri,r 2 ,r 3 ) e are the coordinates of the particle in the inertial reference 
frame, and (e x , e 2 , e 3 ) are the basis vectors fixed in the inertial reference 
frame. Similarly, r p can be expressed in another reference frame as 

r p = J'i^i + r 2^2 + ^363 (2.2.2) 

where (r!,r 2 ,r 3 ) 6 are the coordinates of r p expressed in the b reference 
frame, and (6 X , 63,63) are the basis vectors of an arbitrary moving reference 
frame. 



Fig. 2.1 Position Vector in Three-Dimensional Space 
Since ( 2 . 2 . 1 ) and (2.2.2) describe the same vector, the components of the 
e reference frame must evidently be related to the b reference frame. L he 
relation between the two bases can be established by writing the orthogonal 
projection of r p with respect to the e basis vector so that 

r\ - r p • e x = (e x • 6 x )r x + (e x • 6 2 )r 2 -r (e x • 6 3 )r 3 


( 2 . 2 . 5 ) 
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r 2 — r p “ e 2 — (e 2 • ^ 1)^1 + (e 2 * 6 2 )r 2 + (e 2 * ks ) r 3 (2.2.4) 

r 3 = r p ’ —3 = (^3 ‘ ^l) r l + (^3 * —2)^*2 “I" (—3 ‘ — 3) r 3 (2.2.5) 

These relations can be written in the following matrix form: 


r e = R r 


Th 


(2.2.6) 


where 


R J 


«1 


*1 

—2 

«i 

—3 

—2 

■*1 

—2 

—2 

«2 


.—3 

■*1 


’ —2 

«3 

‘ ^3 


( 2 . 2 . 7 ) 


r e = [rj , r 2 , r!] T , and r 6 = [rj,r 2 ,r 3 ] r . Equation (2.2.6) can be explicitly 
rewritten as the relation of the two basis vectors: 


e = R r b (2.2.8) 

where matrix R is called the coordinate transformation matrix or direction 
cosine matrix between the two sets of axes. 

Since r p preserves the property of constant length regardless of the 
basis vectors selected, 


r„ = r 


er R r Rrf, - 


e T e 

r r 

v p 


which implies that 


R t R = I 


( 2 . 2 . 9 ) 


( 2 . 2 . 10 ) 


or 

R -1 = R r ' (2.2.11) 

A matrix satisfying relation (2.2.11) is called an orthogonal matrix. Pre- 
multiplying (2.2.8) by R and recalling (2.2.10) yields the relations of an 
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arbitrary moving reference frame that is expressed in terms of the inertial 
reference frame as 

b = Re (2.2.12) 

2.3 Angular Orientations of a Particle 

From the previous derivation, we conclude that at any specific time 
the position of a particle, which may be expressed in terms of suitable sets 
of reference frames, can be specified by a transformation matrix. As time 
passes, the position vector orientation changes and so does coordinate trans- 
formation matrix. This leads to the development of the Euler theorem which 
provides us with a foundation to develop various types of angular orien- 
tations so that the coordinate transformation matrix may be determined. 
Euler’s theorem states that two arbitrarily oriented dextral basis vectors b 
and e, with common origin can be made to coincide with one another by 
rotating one of them through a certain angle about an axis which passes 
through that origin. In short, any rotation can be described by rotating a 
vector about a proper unit axis n through an angle 4> as shown in Fig. 2.2. 
The rotational operator acting on the vector can be represented by 

R(n, <f)) -nn + cos 0(1 - nn) + sin<pn x I (2.3.1) 

Note that if the full coordinate transformation matrix is used, it means that 
we choose to parametrize R by nine parameters, the nine direction cosines 
themselves. However, the orthonormal property of the 3 \ 3 coordinate 
transformation matrix R will lead to the consequence that it can only be 
represented by three independent parameters if the six orthonormal con- 
straints are imposed. The choice of these parameters presents an important. 
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aspect in computing the angular orientation of bodies in MBD systems. 
Which kind of parameters one should adopt are a matter of judgment by 
individual researcher. The fewer the parameters, the fewer constraint equa- 
tions need to be satisfied. However, sometimes it is necessary to compute 
the coordinate transformation matrix at every time integration step, which 
cancels some of the advantages by using a small number of parameters. 
Moreover, the use of three parameters always lead to a singular coordinate 
transformation matrix when certain angles are reached. To specify R with 
various parameters, several commonly used parameters will be listed along 
with some important properties so that the kinematic relationships of the 
parameters and their corresponding rotational operator are defined. 



Fig. 2.2 System Rotational Coordinate 
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2.4 Euler Angles 

The most common minimal set of parameters used to describe the 
angular orientation of a body in space are Euler angles. These angles provide 
a set of coordinates without involving any constraint equations. The Euler 
angles formalism consists of three successive rotations, started by rotating 
the A:- ax is through angle a of a specified orthogonal basis \i,j,k}. The 
resulting coordinates are labeled Next, rotate the i '- axis by an 

angle 0 so that another set of coordinates are obtained. Finally, 

rotate the fc^-axis by an angle 7 to produce the desired system rotational 
axes. To express the effective rotational axis n through an angle 4> in terms 
of these three successive rotations, the rotational operator can be written 


R(n ,<f>) = R(fc",7) -R(*',/?) -R(k,a) (2.4.1) 


or 


R(n» 


c^ca — s~ic(3sa s~jca + cqcfisa s0sa 
— c^sct — s'ycficet —s'ysot + c")c0ca s(3ca 
sqs/3 — c'fsfi c0 


The successive rotations in (2.4.1) are 


b" = R (fc",q)b' 
b' = R(i',/J)b 


b = R (k, «)e 


with (c cos , s — sin), 


cq s') 0 

—S' 7 C 7 0 

0 0 1 


( 2 , 1 . 2 ) 


(2,1.3) 


R (k",l) = 
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R(*',0) 


R(fc, a) 


10 0 
0 c/3 s(3 

0 —s(3 c/3 

cot sa 0 
—sa cot 0 
0 0 1 


(2.4.4) 


Multiplying these together yields (2.4.2). The coordinate transformation 
matrix in terms of Euler angles presents some numerical difficulties: first, 
they involve trigonometric functions which are numerically expensive to 
compute; second, the coordinate transformation matrix becomes singular 
when (3 = nir y n — 0, ±1, ±2, ..., in which case both rotations along k and 
k" become collinear. This can be illustrated by setting (3 = 0 so that 


R = 


c(a + 7) s(a + 7) 0 

— s(a + 7) c(a + 7) 0 

0 0 1 


which represents a single rotation a + 7 about A:-axis. 

In mechanical analysis, sometimes it is necessary to calculate Euler 
angles by using a given coordinate transformation matrix where 


R(a,/?,7) = 


r 11 

r 12 

ri3 

^21 

r 22 

^23 

^31 

r 32 

r 33 


To determine the corresponding Euler angles, we calculate first 


(2.4.5) 


a = 



(2.4.6) 


by recalling (2.4.2) so that (3 and 7 can be evaluated without any ambiguity 


s(3 = rizsa. + r 2 3 ca ; cj3 = r 33 


(2.4.7) 


C 1 — r \\ ca - — r 2 i s l i -S7 - r 13 ca - r 2 2 $~i 
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2.5 Rodrigues Parameters 

The Rodrigues parameters are defined as 

7 = n tan — 

' 2 


(2.5.1) 


where n and (p are the unit vectors along the rotational axes and the rotation 
angle. Obviously 


so 


that 


T 4 . 2 & 

* 7 = tan — 
1 1 2 


2 4 > 

COS — 


Since 


2 1+77 


2 ^ 

cos <f> — 2 cos — — 1 


(2.5.2) 


(2.5.3) 


(2.5,1) 


substituting (2.5.3) into (2.5.4) yields 


COS (p = 


1 - 7 r 7 
1 + 7 T 7 

Again, nsin^ can be written in terms of 7 as 


(2.5.5) 


. , _ . <t> <t> 2 <P 2*7 

n sin 0 = 2n sin - cos — = 27 cos - = — 


7 T 7 


(2.5.6) 


Replacement of nsincji) and cost?!) in (2.3.1) by (2.5.6) and (2.5.5) leads to 


R = 


l + l 1 ! 


[(1 - 7 T 7 )I + 2(77 r - 7)] 


(2.5.1 


where 


0 -73 72 

7 = | 73 0 - 7i (2.5.8) 

-72 7i 0 . 

As in the case of Euler angles, the Rodrigues parameters use a minimal 
set of three variables. Unlike the Euler angles, the Rodrigues parameters 
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expression do not involve trigonometric functions. This can be an advantage 
for actual numerical computation. However, they have the disadvantage 

of becoming infinite if the rotation angle (f> becomes ±kn,k — 1,2, 

Again it is sometimes desirable to compute the Rodrigues parameters if the 
coordinate transformation matrix is given. These relations are obtained by 
subtracting symmetrically from the off-diagonal terms of (2.5.7) which yield 

4 m = (1 + tan 2 ^}{Rjk ~ Rkj) (2.5.9) 

Since 

1 + tr(R) = " , , (2.5.10) 

1 + tan ^ 

therefore by substituting (2.5.9) into (2.5.8), we obtain the expression of 
computing 7, that 

7t = 1 + ir(R ) {R]k ~ Rkj) (2 ' 5 ' U) 

Note that if 1 + £r(R) = 0, then n approach infinity which occurs when 
< f> = ±k 7r as is concluded from previous definition. 

2.6 Euler Parameters 

To avoid degeneration of the coordinate transformation matrix for 
certain values of the parametrizing variables, one has to use more than the 
minimal set of three parameters. The choice of the parameters to represent 
the angular orientations of bodies in MBD systems needs to satisfy the 
following requirements: 

(1) Singularity should not occur for any chosen parameters. 

(2) To prevent expensive calculations of trigonometric functions, an alge- 
braic description of finite rotations is preferred. 
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(3) To avoid redundancies in descriptions of parameters, a minimal set of 
parameters is preferred. 

In the present research, Euler parameters have been chosen to represent the 
angular orientations of the bodies for the following reasons: 

(1) Euler parameters satisfy the singularity free property that other sets of 
rotational parametrizations such as Euler angles can not provide. 

(2) Euler parameters preserve the algebraic description of finite rotations 
of bodies in the systems. 

(3) The use of Euler parameters may drastically simplify the mathematical 
formulation. 

Euler parameters are defined as 

<f> 

q 0 = cos - 

, (2-6.1) 
. <P 

q„ = n sin - 

Mn 2 

with the constraint equation 

9o + 9i + 92 + ?3 = 1 (2.6.2) 

where q n = (<7i , < 72 , <73| T - The time derivative of (2.6.2) is given by 


qoqo + = 0 


(2.6.3) 


Introducing the standard trigonometric relationships 

cos 6 = 2 cos 1 

2 

. 6 4 > 

sin 6 = 2 sin — cos — 

Y 2 2 

and substituting (2.6.1) into (2.3.1), the rotational operator dyadic R be- 


( 2 . 6 , 1 ) 


comes 


R = (2 ql - 1)1 + 2(q„q£ - q 0 q n ) 


(2.6.5) 
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or 

R = 2 

where 


‘90 + 01-5 

0102 - 0003 
_ 0103 + 0002 


Qn = 


0102 + 0003 
<7o + 02 2 - 5 
<72 93 — 0001 

0 -93 

93 0 

92 9i 


9i 93 ~ 9o92 
9293 + 9o9i 
<7o + 9 3 2 -|. 

92 

-9x 

0 


( 2 . 6 . 6 ) 


(2.6.7) 


In contrast to Euler angles and Rodrigues parameters, the coordinate trans- 


formation matrix (2.6.6) can not become singular. Again, if R is given, the 
corresponding Euler parameters can be determined by taking the trace of 


R from (2.6.6) so that 



1 + fr(R) 
4 


( 2 . 6 . 8 ) 


Substituting (2.6.8) into the diagonal terms of (2.6.6) results in 



1 + 2 Ra - ir(R) 


i 


1,2,3 


(2.6.9) 


Equations (2.6.8) and (2.6.9) determine the magnitudes of the Euler pa- 
rameters. The off-diagonal terms of (2.6.6) can be used to decide the sign 
of the Euler parameters. Subtracting and adding symmetrically from the 
off-diagonal terms of (2.6.6) yields 


and 


91 

92 

93 


R 32 — R 23 R 32 — R'li 

^ 0q - 

400 4 0 ! 

Rl3 ~ R3 1 R\3 — ^31 

40o 402 

f ?21 — ^12 R -21 ~ R \2 

A ’ 90 = . 

400 4 93 


R 12 + ^21 

0102 — 

4 

R23 + Rz2 

9 2 03 = ; 


( 2 . 6 . 10 ) 


(2.6.11) 
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i ?31 + Rl 3 

g 3 gi - - 

According to previous derivations, the following algorithm is used to deter- 
mined the Euler parameters when R is given [44]: 
det = max ( tr(R), Rn, R 22 > R 33 ) 
if ( det = tr(R) ) then 
use (2.6.8) to find q 0 
use (2.6.10a) to find 91, 92, and 93 
elseif ( det = R t i ) then 
use (2.6.9) to find q t 
Compute q 0 = ^ Rk] ^ lk ' > from (2.6.10b) 

Compute q p = p # i from (2.6.11) 

endif 

2.7 Angular Velocity 

Consider the orientation of the b basis with respect to the e basis 
as given by (2.2.12). The time derivative of b is 

b = Re + Re (2.7.1) 

Since e is a fixed basis vector, which implies e = 0, therefore 

b = Re = RR T b (2.7.2) 

To relate R and R, we differentiate the identity matrix (2.2.10) with respect 
to time: 

R T R + R r R = 0 (2.7.3) 

By assuming that R = SR and substituting into (2.7.3) we get 


R r S r R + R t SR = 0 


(2.7.4) 
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Premultiplying (2.7.4) by R and postmultiplying by R r yields 

S + S r = 0 (2.7.5) 

which implies that S is a skew-symmetric matrix. Hence the matrix kine- 
matic equation for the rotation can be defined as 

R = SR = u) r R (2.7.6) 

or 

Cj t = RR t (2.7.7) 

where 

0 — CJ 3 U>2 

Co = CJ 3 0 — to 1 (2.7.8) 

— U2 (jJ\ 0 

The three components in (2.7.8) are the angular velocity components of the 
moving b basis relative to the inertial e basis that can be written into the 
following form by substituting (2.7.6) into (2.7.2) 

b = RR T b = ii r RR T b = w r b (2.7.9) 

where the angular velocity vector, u can be written as 

00 — Uibi + (^ 2^2 “h ^sb 3 (2.7.10) 

From the present derivation, we conclude that the angular velocity u ; is a 
function of the coordinate transformation matrix R and its time derivatives. 

2.8 Time Derivatives of Euler Parameters 

In this section, the relations between Euler parameters and angular 
velocities are derived. These relations are established by taking the time 
derivative of (2.6.5): 

R = 4 < 7 o< 7 oI + 2 q n q£ + 2q n q£ - 2 < 7 0 q n - 2 < 7 0 q n 


( 2 . 8 . 1 ) 
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Substituting (2.6.5) and (2.8.1) into (2.7.7), the angular velocity, w, can be 
expressed in terms of Euler parameters and their time derivatives as 


where 


W = 2<7 0 q n - 2q 0 q n ~ 2q n q ri = 2Tq 



-01 

<7o 

Q 3 

-Q 2 

T = 

-<72 

<73 

Qo 

Q l 


. “93 

<72 

~Q l 

Qo 


Differentiating (2.8.2) with respect to time yields 


Cj = 2Tq + 2Tq 


(2.8.2) 


(2.8.3) 


(2.8.4) 


Expansion of the product Tq shows that it vanishes, and so does Tq. Hence 
oj = 2Tq and its inverse relation is 

q = -T r w - -( wTw )Q (2.8.5) 

2 4 

Note that the scalar oj t oj = u> 2 can also be written as 4q r q = w 2 if (2.8.2) 
is used. Appending the differential form of the constraint equation (2.6.3) 
to (2.8.2), the angular velocity can be written in terms of Euler parameters 


as 



Qo 

Q l 

Q 2 

93 

( Qo 

~Q l 

Qo 

Q 3 

— 92 

U, 

~Q 2 

~Q 3 

Qo 

9i 

Q2 

_~Q 3 

Q 2 

-Qi 

9o . 

[ Q3 


The inverse of (2.8.6) is 




" 0 

-wi 

— 

— CJ3 

[ 90 1 

1 


0 

CJ 3 

— CJ2 

r 1 

2 

U)2 


0 

^1 

1 Q2 


_^3 

CJ-2 

-U/! 

0 

[ Q3 I 


= d(w)q 


(2.8.6) 


f'2.8 




where q = [qo, qi , < 72 , <?3] T - In chapter 5, these derivations will be used to 
formulate the computational sequences and obtain the angular orientations 


of bodies in the system. 
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2.9 Velocity and Acceleration of a Particle 

In the previous section, we have studied the rate change of a vector 
fixed in the moving reference frame. In the present section, we derive the 
expression for the time derivative of a vector whose components along the 
moving frame are varying with time. Such a vector can be expressed in 
terms of two different basis vectors as shown in ( 2 . 2 . 1 ) and ( 2 . 2 . 2 ) where 

r p = r p = r % + r 2 b 2 + *363 (2.9.1) 

Differentiating (2.9.1) with respect to time and making the use of ( 2 . 7 . 9 ) 
yields 

r e p = + w x r b p (2.9.2) 

where r b p denotes the time rate relative to b basis and u: x r b p denotes the 
time rate of r p due to the rotational motion of the moving frame. Note that 
(2.9.2) represents the time derivative of the position vector r p in an inertial 
reference frame whereas the vector r p is expressed in terms of a moving 
reference frame that is valid for any vector in space. Thus, differentiating 
(2.9.2) with respect to time, we obtain an expression for the acceleration of 
point p : 

r' = fp + 2u x fp + u x Tp + w x w x (2.9.3) 

where r e p is the acceleration of p in the e basis, is the acceleration of p 
in the b basis, 2 u> x is the Coriolis acceleration, oj x v p is the angular 
acceleration in the b basis, and w x u x r* is the centripetal acceleration. 

2.10 Velocity and Acceleration of a Rigid Body 

Having derived the kinematics of a particle in space, it is appro- 
priate to study the velocity and acceleration of a rigid body. Because an 
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unconstrained rigid body possesses of six degrees of freedom, it is generally 
convenient to choose six coordinates that consist of three translations of a 
point within the body and three rotations about that point, to describe the 
motion of the body in space. To this end, consider a position vector r p (Fig. 
2.3) on the rigid body which can be decomposed to 



Fig. 2.3 Translation and Rotation of a Body 
in Three-Dimensional Space 

r p = r 0 + s p = r r e + l r b (2.10.1) 

where r 0 is the position vector from the origin of the inertial reference frame 
to the origin of the body-fixed reference frame, s p is the position vector from 
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the origin of the body-fixed reference frame to a point p of the body, e is 
the basis vector of the inertial reference frame, and b is the basis vectors 
of the body-fixed reference frame that describes the orientation. Also r is 
the position vector of point o in the inertial reference frame, and 1 is the 
position vector of point p in the body-fixed reference frame. By adopting 
(2.9.2) and (2.9.3), the velocity and acceleration vectors of point p can be 
expressed as 

r p = r + s p + w x s p (2.10.2) 

and 

f p = r + 8 P + 2w x s p + w x s p t w x u x s p (2.10.3) 

Since there is no relative motion between the particle at point p and o for 
the rigid body, s p = 0 and s p — 0. The final velocity and acceleration of r p 
can be expressed as 

r p = re + l T b - r T e + l T a> T b (2.10.1) 

? p = r r e + l r a) T b + l r u> r <I> r b (2.10.5) 

Note that if points o and p coincide, which implies 1 = 0, we can derive the 
equations of motion by separating the translational and rotational motion 
so that different numerical algorithms may be applied accordingly. 

2.11 Concluding Remarks 

Spatial kinematics relations needed to calculate quantities such as 
the position, velocity and acceleration vectors of particles and bodies have 
been reviewed. In discussing the motion of a particle as being attached to 
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a rigid body, a frame of reference must be specified so that the dynamical 
equations of the body can be derived. In the present derivation, an inertial 
reference frame is used to locate the position of the center of mass of a body. 
A body-fixed reference frame is then employed to locate the position of a 
particle in the body. Such hybrid reference frames are chosen to decouple 
the translational and rotational equations so that an efficient numerical 
algorithm can be formulated as discussed in chapter 5. 

Three representations of rotation have been studied. The advan- 
tages and disadvantages of these angular representations are discussed. Eu- 
ler parameters have been chosen in the present derivation because they lead 
to simple algebraic equations that do not require the evaluation of trigono- 
metric functions and they give a singularity-free representation of rotations. 
Other angular representations may require trigonometric functions and/or 
suffer from singularity for certain parameter values. 



CHAPTER III 


EQUATIONS OF MOTION FOR MBD SYSTEMS 
3.1 Introduction 

This chapter deals with the formulation of the equations of motion 
for MBD systems using variational methods. There are several advantages 
of employing variational methods in dynamics. First, the system of particles 
and rigid bodies is considered as a whole rather than being separated into 
individual components. Second, dynamic problems are formulated in terms 
of kinetic energy and work, both of which are scalar quantities. Third, 
constraint forces do no work. Fourth, the use of generalized coordinates 
makes the formulation versatile. In this regard, the reference frames and 
velocities reviewed in chapter 2 will be adopted to describe the configuration 
and motion of bodies in MBD systems. 

As indicated in chapter 2, an inertial frame is used to described 
the translational motion whereas a body-fixed frame is used to described 
the rotational motion of bodies in the system. This frame decomposition 
causes the mass matrix to be decoupled into translational and rotational 
equations. This kinematic representation is introduced into t he principle of 
virtual work to obtain dynamic relationships between the constraint forces 
and their kinematic constraint conditions and thus produce the governing 
equations of motion. When d’Alembert’s principle is used in conjunction 
with the principle of virtual work, we extend the principle of virtual work 
to dynamic systems that are composed of an arbitrary number of rigid bod- 
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ies with an arbitrary number of constraints that are used to restrict the 
motion of the bodies. The present formalism considers the motions of indi- 
vidual bodies as initially independent, and then applies restriction on those 
motion by the introduction of kinematic constraints. Such constraints are 
incorporated through the method of Lagrange multipliers. The resulting 
system of equations, which consists of second-order differential equations 
that introduce Lagrange multipliers as constraint forces as well as the alge- 
braic constraint equations as constraint conditions, are known as differential- 
algebraic equations (DAEs). To cover further development in flexible MBD 
systems, the equations of motion that include elastic deformations, which 
have been derived by Downer [52], are given in time discrete form by taking 
the advantage of the previously chosen reference frames and formulation. 

3.2 The Principle of Virtual Work 

Since an MBD system involves a number of interconnected bod- 
ies, the study of its dynamics is simplified in many respects by considering 
the system as a whole rather than as a collection of components obeying 
Newton’s laws of motion. This is accomplished, as noted previously, by 
basing the derivation on an overall scalar quantity: generalized work. Con- 
sequently, the principle of virtual work will be used to establish the system s 
equilibrium conditions. This principle may be stated as follows: 7 he work 
done by all the forces acting on a system in static equilibrium, during a 
virtual displacement compatible with the constraints of the system is equal 
to zero. The mathematical expression is 

n 

sw = • 6n = 0 (3.2.1) 

1=1 

where n denotes the total number of bodies, <5W denotes t ho virtual work 
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of the system, P t denotes the resultant forces acting on each body, and 6t{ 
denotes the virtual displacement of each body. 

To interconnect and restrict the motion of bodies in the overall 
system, kinematic constraints are imposed. Two types of constraints must 
be distinguished: 

(1) Holonomic: constraints that depend only on position. Such kinematic 
restrictions may be expressed as algebraic relations: 

<S> h {r p ,t)=0 (3.2.2) 

The variation of the holonomic constraints is given by 

d&h 

6& h = — =<5r p = B h 6r p = 0 (3.2.3) 

OT p 

(2) Nonholonomic: constraints that depend both on position and velocities. 
Such kinematic restrictions are expressed as in differential relations: 

^n/i(rp } rp,£) = Bix/jTp = 0 (3.2.4) 

The variation of the nonholonomic constraints is given by 

S®nh = B n /^rp = 0 (3.2.5) 

where h and nh refer to holonomic and nonholonomic constraints. Hence 
when systems are subjected to constraints, one may separate the resultant 
forces P t into applied forces F^ and constraint forces F£, so that 

P: = F* + F c t (3.2.6) 

Substituting (3.2.6) into (3.2.1) yields 

n n 

6W = • 6n + • 6n = 0 (3.2.7) 

1=1 1=1 
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Since the variations <5r t do not violate the prescribed constraint forces, the 
work performed by the constraint forces in any virtual displacement is equal 
to zero, therefore we conclude that 

n 

■ Sr i = 0 (3.2.8) 

»=i 

Note that for systems with constraints, the virtual displacements <5r, are 
not all independent. Thus we cannot interpret (3.2.8) as F^ = 0. For an 
unconstrained system, the principle of virtual work can be used to calculate 
the equilibrium position of the system as 

n , n . dV 

<5W = F“ ' Sr t = SV = £(— • 6 r t ) = 0 (3.2.9) 

i=l i=i 1 

where V is the potential energy. Since by hypothesis the virtual displace- 
ments <5r, are all independent the static equilibrium conditions can be ob- 
tained as expected: 

F“ = = 0 (3.2.10) 

dri 

If a system is subjected to holonomic constraints 

*(r)=0 (3.2.11) 

the method of Lagrange multipliers is used to augment the potential energy. 
According to this method, we multiply each of the constraints (3.2.11) by 
an undetermined multiplier Ay, and add all resulting expressions to the 
potential energy V to get 

m 

V a = V + £(A,- • $j) (3.2.12) 

/= i 

where V a is the augmented potential energy and m denotes the total number 
of the constraint equations. The variation of the augmented potential energy 
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subjected to the constraint condition (3.2.11) can be written 

n Qy m 

*V = Efc ■ «■-.•) + ' ‘*i) = 0 (3-2.13) 

i=i * y=i 

Substituting (3.2.3) into (3.2.13) yields 


,3V v— v d$j. . 

(^T + • -^f) • <5r « = °. 1 = 1, ••■>» (3.2.14) 

* y=i 1 

Note that the virtual displacements Sr, in (3.2.14) are still not independent, 
but the m values of Ay can be chosen so that 


3V d&j 

— — +> Ay — — = 0, i = n — m + l,n — m + 2, ...,n 

3r, ' 3r, 

;=i 


whereas the remaining n — m virtual displacements Sr,(t = 1, .. 


be treated as independent variables so that 


(3.2.15) 
n — m) can 


3V JV 3$, 

ss + E^-fcf" 0, , = 1 "" 

J = 1 


n — m 


(3.2.16) 

From (3.2.15) and (3.2.16), we obtain the following equilibrium conditions 


■ -sf = o. * =1 » < 3 - 217 ) 

J=1 

This procedure enables us to treat all the virtual displacements as indepen- 
dent variables by expanding from n — m unknowns to n + rn unknowns with 
n values of r x and m values of Ay. Now, recalling (3.2.10) and comparing 
the expression of (3.2.7) and (3.2.17), we arrive at the conclusion that the 
system equilibrium conditions are enforced by the presence of the constraint 
forces 


£v 

j = i 


3$ 

3r, 


i = 1 n 


( 3 . 2 . 18 ) 
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This important result provides the relationship between the constraint forces 
and the kinematic conditions of the bodies in the system. The principle 
of virtual work was originally stated for a system in static equilibrium. 
Nevertheless, the principle can be applied to dynamic systems by simple 
recourse to d’Alembert’s principle, which gives the dynamic equilibrium by 
including the inertial forces of the system with constraints. 

3.3 D’Alembert’s Principle 

D’Alembert’s principle states that the law of state equilibrium ap- 
plies to a dynamic system if the inertial forces as well as the external and 
constraint forces are considered as applied forces acting on the system. Thus, 
for a body with density p, the dynamic equilibrium condition is given by 

jp‘_F°_F c = 0 (3.3.1) 

where F* = pr p are the inertia forces, and r p are the acceleration vectors. If 
we apply d’Alembert’s principle in conjunction with the principle of virtual 
work, the principle of virtual work is extended by writting the following 
equation: 



Where F a may be considered to include many types of force acting on 
the body: viscous forces which resist velocities; spring forces which restore 
position equilibrium; and independent defined external forces. We shall refer 
to (3.3.2), which includes both the principle of virtual work and d’Alembert’s 
principle, as d’Alembert’s principle of virtual work. In present chapter, we 
use this formulation to derive the equations of motion for MBD systems. 
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3.4 Governing Equations of Motion 

To derive the equations of motion for MBD systems, we start with 
an unconstrained rigid body by using (3.3.2) where F c = 0. There are 
many possibilities of choosing 8 r p depending upon the coordinates one has 
employed. In present derivations, we adopt the velocity and acceleration 
vectors of point p derived in (2.10.4), and (2.10.5). The virtual displacement 
6r p can be obtained as 

6r p = 6r T e + l T «5b = 6r T e + l T <5a T b (3.4.1) 

and the virtual rotational tensor 6a is 


6a = -6 RR J 


(3-4.2) 


Substituting these two equations into (3.3.2) yields 

f {6r T e + 6a T lb) ■ [p(? T e + l r <i> r b + l T w T w T b) - i T b\dV = 

Jv 

<5r r [M(f 4- R r r Jcj -f R r u7r Ju>) — F] + <5a: r [M?J Rr + Jc5 + u3u - M 0 J = 0 


(3.4.3) 


where 


M 


= f pdV, 3 = f 
J v J v 


pdV, J= / p\Y dV, Mf c = f pldV 
l V J V J V 


F = / 

J V 


R f d\\ M 0 
v Jv 


-L 


(3,1,1) 


MdV 


Performing the variation of <5 r and 6a independently, the equations of mo- 
tion for a unconstrained rigid body can be written in the following forms 


(3.4.5) 


M(f + R T rJd> + R r u>rju>) - F = 0 
Mr ^Rr + Jc5 + CiJu — M 0 = 0 


(3.4.6) 
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A considerable simplification can be made in the equations of motion if the 
body-fixed coordinates are chosen such that the principle axes coincide with 
the center of mass. With this choice, all products of inertia vanish since 

r c = 0 (3.4.7) 

and (3.4.5) and (3.4.6) reduce to 

Mr = F (3.4.8) 


J Cj + ujJui - M ( J 


(3.4.9) 


Note that the translational and rotational mass matrices can be expressed 


as follows: 


M = 


m 

0 

0 


0 

m 

0 


0 

0 

m 


(3.4.10) 


J = 




'2,3 


(3.4.11 


*A,1 j 1,2 
^ 2,1 ^ 2,2 
L>^3,1 J3,2 Jz,3 

where M and J denote the mass and moment of inertia of the body. Equa- 


tions (3.4.8) and (3.4.9) are known as Newton-Euler’s equations of motion. 
Euler equations (3.4.9) are widely used in solving for the rotational motion 
of a rigid body. Note that, however, they are in general nonlinear and it 
may be difficult to solve analytically for angular velocity u; as a function of 
time. Furthermore, the time integral of u does not correspond to any phys- 
ical rotational representation that can be used to describe the orientation 
of the body. So if one wishes to find the angular orientation of a body, a 
set of parameters must be chosen in order to find the relation between the 
parameters that orient the bodies and their corresponding angular velocity. 
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For MBD systems, the constraint conditions are introduced into 
d’Alembert’s principle of virtual work via Lagrange multipliers to restrict 
the motion of the bodies in the system. From a formulation point of view, 
there are several ways to impose the kinematic relationships between bodies 
during the motion. In the present derivation, we use the description of the 
unconstrained motion to describe each of the bodies separately. Therefore, 
the virtual displacement of (3.4.1) is not a kinematically admissible one for 
the constrained systems. The method of Lagrange multipliers must then 
be introduced to incorporate the constraint conditions into d’Alembert’s 
principle of virtual work as has been indicated in the previous sections. To 
apply this method the constraints are multiplied by undetermined Lagrange 
multipliers A and added to the virtual work of the unconstrained system: 



f) + $* • A \dV = 0 


(3.4.12) 


where 6r p , p,f p ,f , and dV are defined in the previous derivations, A are the 
Lagrange multipliers and <5$ are the variations of the constraint equations. 
The augmented terms represent the work of the constraint forces, provid- 
ing the reaction forces which are exerted on account of given kinematica! 
constraints. 

Substituting (3.2.3) and (3.2.5) into (3.4.12) yields 


/ [<5r p • (pf p - f) + 6® h ■ A h + ■ A nfe jdV' 

Jv 

f fop ■ {pip -f + Bl\ h + B* h \ nh )dV = 0 


(3.1.13) 


Performing the variation of 6r and 6 a independently, the equations of mo- 
tion for MBD systems can be derived from (3.4.13) in the following matrix 
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form 


M 0 
0 J 




■ f + 


.{ p . \ 

\ M c — wJo; J 


(3.4.14) 


Augmenting (3.4.14) with the constraint equations (3.2.2) and (3.2.4), the 
differential-algebraic equations result: 


Mu + B T A = F 


(3.4.15) 


that are subjected to satisfy holonomic constraints 


®fc(u,i) = 0 


(3.4.16) 


and nonholonomic constraints 


6 n ^(u,u,i) = BnfcU = 0 


(3.4.17) 


where ii = [r,w] r , B is the gradient of the holonomic and nonholonomic 
constraints (or constraint Jacobian matrix), A is its corresponding constraint 
forces, F is the forces that include external forces and inertia forces due to 
centrifugal acceleration, and u is the generalized displacement vector. The 
mass matrix for j-th body is given by combining (3.4.10) and (3.4.11) as 


M J = 


m 

0 

0 

0 

0 

0 

0 

m 

0 

0 

0 

0 

0 

0 

m 

0 

0 

0 

0 

0 

0 

Ji.i 

J l,2 

Jl, 3 

0 

0 

0 

J2.1 

j 2,2 

*^2,3 

0 

0 

0 

J 3,l 

Jz,2 

*^3,3 . 


(3.4.18) 


and the force vector for j-th body is 


F : = 


/_ 

M 0 - uiJlj 


(3.4.19) 
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In the present derivation we have replaced a system having n — m unknowns 
by one with n + m unknowns, considering the Lagrange multipliers A as 
additional variables where n is the total number of degrees of freedom be- 
fore imposing constraints and m is the number of constraint equations. The 
advantages of the present derivation are: first, the mass matrix as shown 
in (3.4.18) can be partitioned into translational and rotational sets of equa- 
tions, which later will lead to a convenient computational algorithm treats 
the rotational equations and the translational equations with different proce- 
dures. Second, the method of Lagrange multipliers preserves the symmetry 
of the resulting equations for all coordinates without distinguishing between 
dependent and independent variables. Third, the constraint Jacobian ma- 
trix that defines the kinematic relationships between interconnected bodies 
can be generated by using a set of stand-alone joint modules. Fourth, the 
presence of closed loops in the system topologies, require no special treat- 
ment so that preprocessing to identify independent variables can be avoided. 

The velocity and acceleration equations for holonomic constraints 
are given by 


6^ = B^u + & t (3.4.20) 

= B h u + B h u + 2& ut u + & tt (3.4.21) 

The acceleration equation for nonholonomic constraints is given by 

= B n/l ii + B n fcU + 2$ u tii + $ t t (3.4.22) 

Regardless of the nature of the constraints, the equations of motion with 
the constraint acceleration equation can be augmented into the following 
matrix form: 


[» T] {:}-('} 


(3.4.23) 
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where c = -(Bu + 2$ u tU + &tt)- Since the left hand side of (3.4.23) 
is symmetric and sparse, several research groups have developed solution 
procedures tailored to solve these constraint-augmented equations. These 
solution procedures will be discussed in chapter 5. 

3.5 Interaction Equations for Rigid and Flexible Bodies 

Up to now, the bodies that comprise the MBD system have been 
rigid. This assumption does not hold when the bodies in the mechanical 
system are subject to elastic deformation that must be taken into account. 
The formulation presented in this section has been motivated by further 
developments in analysis and design of large-scale systems that consist of 
interconnected rigid and flexible bodies, all of which may undergo large an- 
gular rotations as well as deformation. As discussed in previous sections, the 
bodies (rigid or flexible) in MBD systems are treated initially independent of 
each other. Kinematic relationships between adjacent rigid or flexible bod- 
ies are specified through a set of nonlinear algebraic constraint equations 
that depend on the position and time. 

The purpose of this section is to impose these kinematic relation- 
ships into rigid or flexible bodies of multibody systems so that their equa- 
tions of motion, which can be written in time discrete form, can be obtained. 
Since that the formulation of flexible body dynamics is well documented, 
e.g., in Downer et al. [52], only the body interfacing requirements will be 
outlined. There are essentially two different connection cases to be consid- 
ered in flexible MBD systems: first, two flexible bodies are connected by 
a specific joint; second, a flexible body connects a rigid body with a given 
kinematic relationship. These approaches are illustrated in the following 
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sections as initial development for the two-stage staggered explicit-implicit 
algorithm, discussed in chapter 5, which is used to numerically integrate 
these sets of nonlinear equations. 

3.5.1 The Equations of Motion: Interaction of Flexible Bodies 

The discrete equations of motion for this approach can be expressed 
(Downer, Park, and Chiou [52]) as illustrated in Fig. 3.1 where 



Y(e 2 ) 


/ 

Z(e 3 ) 


*■ X(«j) 


Fig. 3.1 Interaction of Flexible Bodies 
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Mii + D(u) + S(u) + B^A h + B^Anh = F 


( 3 . 5 . 11 ) 


or 


Mii + b£ A/i + B^A n /i — Q 

subject to the constraint equations, 

*fc(™,0=° ; ^nfc(u,u,0 = B nh il = 0 

where 


M = 


Mi 0 0 

0 My 0 

0 0 M fc J 


and A 'frib — ..., 

u = [U(*,1), ...,U(t,nt)> W(y il ),...,U(y ) nj),«(fc,l ) 1 ---> * ( fc.ixfc) ] ' 


B h = 


B 

B 

y.i) 

0 

0 

- 

0 

0 

B(j,nj) B(k, 1) 


' F{i,i) 

-Su 

,1) - 

%,D ' 



- s [it 

ni~) 

B(i,n i) 


F U, i) 

- S ti, i) ~~ 

D UA) 

= < 


- 


> 



1 

c 

'—1 

1 



F (M) 

- S { kA) - 

D (k,l) 


i F (k,nk) 

- S(k 

nk) 

~ D(k,nk) . 


( 3 . 5 . 12 ) 


( 3 . 5 . 13 ) 


( 3 . 5 . 14 ) 


( 3 . 5 . 15 ) 

( 3 . 5 . 16 ) 


( 3 . 5 . 17 ) 


In the above equations m, nj and nk are the total number of discrete nodal 
points, subscript ( nb,nd ) denotes the nd - th node of the nb- 1 h flexible body. 
M is the mass matrix of i-th, y'-th, and /c-th bodies, diag are the diagonal 
block matrices of each individual body, D(-) is the generalized velocity- 
dependent force, S(-) is the internal force operator due to member flexibility, 
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is the gradient of the holonomic constraints that connect the nodes with 
prescribed joints, B n ^ is the nonholonomic constraint Jacobian matrix, 
are the holonomic constraint forces, A n ^ are the nonholonomic constraint 
forces, F are external forces, and u is the generalized displacement vector. 

In the present time discrete form, the flexible bodies are initially 
treated as independent of each other. Their kinematic relationships are 
then imposed by given specific constraint conditions at certain nodal points. 
Thus, the Lagrange multipliers will only be computed via the quantities of 
these constraint nodal points. We further address this issue in chapter 5 
where the two-stage staggered explicit-implicit algorithm is developed. 

3.5.2 The Equations of Motion: Interaction of Flexible and Rigid Bodies 
The major difference introduced by the presence of rigid bodies per- 
tains to the construction of the constraint Jacobian matrix, which signifi- 
cantly affects the computation of the constraint forces. The discrete equa- 
tions of motion for this case can be expressed as shown in Fig. 3.2 where 

Mu + D(ii) + S(u) + B^A* + Bl h X nh = F (3.5.21) 

or 

Mii + + B^Ari/i = Q (3.5.22) 

subject to the constraint equations, 

$fc(iM)=° ! $n*(u,u,t) = B n /,u = 0 (3.5.23) 


where 

M, 0 0 

0 Mj 0 
0 0 M k 


M = 


(3.5.24) 
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Z<e 3 ) 


Fig. 3.2 Interaction of Flexible and Rigid Bodies 
or 

r Mu do o o ooo' 

0 ... 0 0 0 0 0 

0 0 M (l>0 0 0 0 0 

M = 0 0 0 My 0 0 0 (3.5.25) 

0 0 0 0 M (M ) 0 0 

0 0 0 0 0 ... o 

0 0 0 0 0 0 m ( m) _ 
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U •••■> U(i,ni)> Uj , ..., U(k,nk)\ 


B = 


Q= < 


B u<nl) 0 
0 B {ji nr) B {ktl) 

\ F (i, 1) ~ _ D {iA) ) 

Fj 

F(k,i) - S {kt i) - D (ktl) 


F(k,nk) B(k,nk) ^(fc.nt) 


(3.5.26) 

(3.5.27) 


(3.5.28) 


In these equations, nt and nk are the total number of discrete nodal points, 
subscript (a, 6) denotes the fr-th node of the a-th flexible body, subscript 
j denotes the j-th connected rigid body, M consists of the mass matrix of 
flexible body i, k and rigid body j, D (•) is the generalized velocity-dependent 
force, S(-) is the internal force operator due to member flexibility, is the 
gradient of the holonomic constraints that connects the m'-th node of t ~ th 
flexible body to the left hand side of the j-th rigid body and the nk- th node 
of A>th flexible body to the right hand side of the j-th rigid body, B n ^ is the 
nonholonomic constraint Jacobian matrix, A^ are the holonomic constraint 
forces, X n ^ are the nonholonomic constraint forces, F( a ^) are the external 
forces, Fk includes inertia forces due to centrifugal acceleration and external 
loads, and u is the generalized displacement vector. 


3.6 Concluding Remarks 

Two methodologies for deriving the equations of motion for systems 
with a number of bodies subject to kinematic constraint conditions may be 
distinguished. The first method makes explicit use of constraint conditions 
so that system dependent and independent variables can be identified and 
eliminated thus reducing the system equations to the number of indepen- 
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dent variables. The second method introduces Lagrange multipliers in the 
equations of motion so that DAEs are obtained. The present research makes 
use of the latter method to generate the system dynamic equations in DAE 
form. Advantages gained by this choice are as follows. First as shown in 
(3.4.21), the symmetry of the DAEs for all coordinates is preserved which 
avoids having to distinguish dependent and independent system variables. 
Second, the constraint Jacobian matrix that establishes the kinematic rela- 
tionships of contiguous bodies can be generated by using a set of stand-alone 
joint modules as presented in the next chapter. Third, the system topology 
whether open or closed-loop, require no special treatment, viz., preprocess- 
ing for a-priori identification of independent variables can be avoided. 

The incorporation of flexible bodies, such as beams, in multibody 
systems is also discussed in this chapter. The computational issues regarding 
these discrete forms will be addressed during the development of the two- 
stage staggered explicit-implicit algorithm. 



CHAPTER IV 


KINEMATIC JOINTS AND FORCE ELEMENTS 
4.1 Introduction 

The equations of motion of MBD systems incorporating elastic de- 
formations have been derived in the previous chapter. It is emphasized that 
the individual bodies are originally treated as independent of each other. 
Kinematic relations that link those bodies are established using constraint 
equations. In the previous chapter, however, the physical interpretation of 
the constraint Jacobian matrix has not been clearly defined. To complete 
the derivation of the equations of motion, the constraint Jacobian matrix 
pertaining to specific mechanical systems must be derived in detail to facil- 
itate the development of a modular computer program. 

A common feature in the construction of these constraint conditions 
is the use of joints to describe the interaction of contiguous bodies in MBD 
systems. Joints may range from rigid connectors, which allows no relative 
motion between two bodies, to devices that allows the relative separation 
of the bodies. Consequently, joint descriptions may involve from zero to 
six degrees of freedom. Two types of kinematic constraints, configuration- 
dependent (holonomic) constraints and velocity-dependent (nonholonomic) 
constraints, are used to describe joints. Spherical, universal, revolute, and 
translational joints provide examples of holonomic constraints. A rigid joint* 
provides an example of nonholonomic constraint. In this chapter, the con- 
straint equations pertaining to a spherical joint, universal joint, revolute 
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joint, cylindrical joint, and a rigid joint are derived. As for other joints such 
as translational and skew joints, the kinematic principles discussed in this 
chapter can still be applied accordingly. 

After the kinematic joints are derived, the force elements such as 
gravity, external forces, moments, actuators, dampers, and springs will be 
incorporated into DAEs in order to complete the assembly of a general- 
purpose computer program. Force elements are discussed after the treat- 
ment of mechanical joints because some of the constraint conditions used in 
kinematic joints can be applied to the force elements thus avoiding unnec- 
essary derivations. 

4.2 Spherical Joint 

A spherical joint is characterized by imposing the equality of the 
positions of two connected bodies at a specific common location. This joint 
allows three relative rotational degrees of freedom during dynamic motion. 
Fig. 4.1 shows two adjacent bodies i and j connected by a spherical joint. 
The center of the spherical joint, called p, can be represented by the body- 
fixed coordinates (6\,6‘ 2 ,6^), and (&{ , 6^) respectively. To restrict the 

relative motion of bodies i and j, the algebraic constraint equations are 
expressed as 

= r, + Sp - Tj - sj, = 0 (4-2.1) 

This relation can also be obtained by applying (2.10.1) as 


e T &. 


T 

r. e 


+ lp,b. 


T 

r e 


-Kfr 


(4.2.2) 


Since is not function of time, if we differentiate <&., once with respect to 
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time, the following relationship is established: 

6 a = B a ii 


(4.2.3) 


or 

e T $3=e T B 5 u (4.2.4) 

where B 3 is the constraint Jacobian matrix and u is the velocity vector 
containing the translational and rotational components of bodies i and j, 
namely 

u = [ r, , W, , r, , w, } T (4.2.5) 



Fig. 4.1 A Spherical Joint 

Equation (4.2.3) shows that if one wishes to obtain B a , we need to extract 
the velocity vector u from the time derivative of the constraint equations 
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and treat the remaining terms as the constraint Jacobian matrix. This can 
be demonstrated by the following algebraic calculations 


i • /p m » *71 m • 

= r t e + lp t b i - fy e - lp ; b ; 


(4.2.6) 


By substituting (2.7.2) into (4.2.6) and using the expression of (2.7.7), we 
obtain 


• T *. = (*?' + ,> 


iT ~.T 


T -T*. 


e (f j + R t - W t lp, fy Ry <2)ylpy) 
The cross product of two vectors a and b is given by 


(4.2.7) 


a x 


b = ab = — ba = — b x a 


(4.2.8) 


Making the use of (4.2.8), (4.2.7) can be transformed to 
e T ** = e T (r,- + - fy - Rfl^y) 

= e T [I, (iJ i R t ) r , —I, (lpyRy) 7 '] 

Comparing the results of (4.2.4) and (4.2.9), we obtain the expression of the 
constraint Jacobian matrix B 3 for a spherical joint where 


r t 

U>i 

U, 


(4.2.9) 


B s = [ I , (1 p ,R,) T , -I , -(lpyRy) 7 * 1 (4.2.10) 

and I denotes the (3 x 3) identity matrix. Similarly, the second time deriva- 
tive of & 3 is given by 

& 3 - B 3 ii + B 3 u 

• T~T . T~T 

= B 3 u + R t lp t w t - Rj \ pj u 3 

= B 3 ii + R t T Uil pi Ui - R JwylpyWy (4.2.1 1 ) 

ft " 

Wi 


= B 3 U + [O.RfwJ^O.Ry’wyipy 
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as 

B, = [ 0 , -(i pi i>iR,) r , 0 , (i„jU j R i ) T 1 (4.2.12) 


and 


u = [ r ,• , u>, 



(4.2.13) 


4.3 Universal Joint 

The universal joint fixes two bodies at an arbitrary position in space 
and allows two relative rotational degrees of freedom during relative motion. 
Fig. 4.2 shows two bodies connected by a universal joint. The constraint 
equations for the universal joint can be expressed as if there were a spherical 
joint between connected bodies with two vectors s, and s ; that are perpen- 
dicular. If two vectors remain perpendicular at all times their dot product 
vanishes. This kinematic relationship can be expressed as 


^ unu — 



Since 


s = s T e = l T b = / r Re 


(4.3.1) 


(4.3.2) 


we conclude that 

s T = l T R ; s = R r /_ 

Replacement of s and s T in (4.3.1.b) by (4.3.3) lead to 


= / t r R,R;/ ; 


Time differentiation of (4.3.4) yields 

4>u - ijRiRjlj + /JR.R^ 


(4.3.3) 


(4.3.4) 


(4.3.5) 



55 



If (2.7.6) is applied, (4.3.5) becomes 


$ u = ifufRiRfLj + ijRiRjujlj (4.3.6) 

Since $ u is a scalar, the transpose of the first term of (4.3.6) gives the same 
magnitude as 

$ u = lj RyR f&ili + LjRiRjujlj (4.3.7) 

Applying (4.2.8), we obtain 

4>u = /jRyRf/fwt + ijRiRj^Uj 

= [ o , Lj RjRfL , o , ljR t Rjl j 
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and the constraint Jacobian matrix B u can be expressed as 


B„ = [ 0 , if RjRf/f , 0 , ZfRiRjlJ 


(4.3.9) 


Following the same procedure as in (4.2.11), the B u is found to be 


B u = [0, if wJ R y R f/f +ij RjRT*£, 0 flj + lj*i Rj*jl 




(4.3.10) 


Hence the gradient matrix of the constraint equations for the universal joint 
can be written as 


B 


unv 


b 3 

B u 


(4.3.11) 


4.4 Revolute Joint 

A revolute joint attaches two bodies in space and allow one rota- 
tional degree of freedom during actual motion. Fig. 4.3 shows two bodies 
connected by a revolute joint. The constraint equations for the revolute 
joint can be expressed as if there were a spherical joint connected two bod- 
ies with two vectors s t and Sy that are always remained parallel to each 
other. Mathematically, their cross product is equal to zero. The constraint 
equations for the revolute joint can be expressed as 


4 > 

^ rev 


$ r „ = S, X Sj 


= 0 


(4,1.1) 


Equation (4.3.2) has concluded that s = R r (, time differentiation of s yields 


s = R T l_ = R T ul = -R T /uj 


(4,1.2) 



57 









58 


to be 

B r = [ 0 , SyR fli , 0 , “S iRjlj ] ( 4 - 4 - 5 ) 

The time derivative of (4.4.5) is found to be 


= [ 0 , SjRfI, + SyR^u)i/j , 0 , S^Ry Ij 


(4.4.6) 


Hence the constraint Jacobian matrix for the revolute joint can be written 
as 

b "” = [b‘] (4 ' 47) 

Notice that if two vectors are to remain parallel at all times, only two con- 
straint equations are needed. Equation (4.4.5) yields three algebraic equa- 
tions, of which only two equations are independent, i.e., one of the equations 
can be derived as linear combination of the remaining equations. A tech- 
nique for selecting the proper set of equations from the overdetermined set 
is to compare the absolute values of each row equation of the gradient ma- 
trix of the constraint equations, and select the two equations that have the 
largest terms. 

An alternative approach to modeling a revolute joint is to set up 
two proper vectors that are capable of representing their kinematical re- 
lationships as a revolute joint. This approach is followed below. Let 
b‘ = [ 6 i, 6 2 , M’ and b J = [ 61 , 62 , 63 ]*' be two triad of orthogonal unit vectors 
attached to bodies t and; respectively (Fig. 4.4). The kinematic constraints 
for this revolute joint can be expressed as 




$rl = 62 • b\ 
$r2 = 6^-6* 


(4.4.8) 
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For any configuration, the constraint Jacobian matrix B r of the revolute 
joint is derived and given by 


4 (R‘) T R'V, 1 
4(r , ) t rv 1 / 


where c{, c' 2 , and C 3 are the components of b* and b- 7 . 


(4.4.9) 



Fig. 4.4 A Modified Revolute Joint 


4.5 Cylindrical Joint 

A cylindrical joint provides one translational and rotational degree 
of freedom to two connected bodies. Fig. 4.5 shows the constraint condi- 
tions for a cylindrical joint. The constraint equations are derived from the 
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condition that vector u, must remain parallel to vectors uy and d: 


® cyi — 


$cl = Ui X Uj = UiUj 1 _ 0 

$ C 2 = Ui x d - Uid J 


(4.5.1) 


where u t - and uy are given directional unit vectors based on their body-fixed 
frames so that the two bodies will slide according to that axis, and d is 
obtained from d = rj — r_j. 



Fig. 4.5 A Cylindrical Joint 

Since (4. 5.1. a) has the same form as in (4.4. l.b), B cl and B cl can be found 
in the same way as (4.4.5) and (4.4.6) . If the first and second time derivatives 
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of (4.5.1.b) are taken, B C 2 and B C 2 can be obtained as 


B c2 = [ u, , dR Jit , u,' , 0 


B C 2 = [ u» , dR t -7 t + dR Juili , -d t , 0 


The final constraint equations for the cylindrical joint are: 


*cyl 


B cl 

B C 2 


(4.5.2) 

(4.5.3) 


(4.5.4) 


4.6 Rigid Joint 

A rigid joint by definition allows no relative motion between two 
bodies. Thus us a nonholonomic constraint that can be imposed as a spher- 
ical joint with no relative velocities among the connected bodies. The above 
statement can be expressed mathematically as following equations 

® rigid = | * j-. . 1—0 (4.6.1) 

l *n = Ui - u, = B rj -u J v 1 

The constraint Jacobian matrix of <& rj is given by 

B r y = [ 0 , I , 0 , -I j (4.6.2) 

Hence the gradient matrix of the constraint equations for the rigid joint can 
be written as 

(4.6.3) 

4.7 Force Elements in MBD Systems 

In section 3.3, different types of forces acting on the bodies have been 
discussed. Forces that are commonly encountered in mechanical systems 
include gravitational forces, actuator forces, damping forces, spring forces, 


B 


rigid — 


B, 

Br, 
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and external forces. In the present section, these forces will be formulated 
and incorporated into DAEs as 

F = F W + F !7 + F a +F (i + F*+F / (4.71) 

where F w , F s , F a , Fj, Fjt and F / denote centripetal, gravity, actuator, 
damping, spring and external forces, respectively. These force elements are 
analyzed in further detail below. 

4.7.1 Gravitational Force 

Since the acceleration of gravity is measured with respect to an 
inertial reference frame fixed in the earth, the gravitational force of a body 
with mass m g can be calculated by the equation 

f g = rn g g (4. 7.1.1) 

where f g is the force created by the gravity and g is the acceleration of grav- 
ity. If we choose a gravitational field acting on the negative z direction of the 
inertial reference frame, the force F 9 that contributes by this gravitational 
field on body i is 

F& = [0,0,— / 3 (,) , 0,0,0] T (4.7. 1.2) 


4.7.2 External Forces and Moments 

Consider a force acting on body i at point p as shown in Fig. 
4.6. The moment of about the origin of the body is 

M 0 (t) - s (l) x /W (4.7.2. l) 

where — l T b ^ is the position vector of point p in the ?-th body-fixed 
reference frame. The contribution of and Mq ' to the force vector F / 
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on body * is 

Ff = [/<*'), Af^J r (4. 7.2. 2) 

If a pure moment acts on body i, then the force vector F/ on body i 
becomes 

F^. 0 = [0,m(* ) ] T (4.7. 2.3) 



Fig. 4.6 A Point Force Acting on a Body 
4.7.3 Actuator Forces 

An actuator is a force element that provides a constant or a time- 
dependent pair of forces acting on two bodies in MBD systems. The direc- 
tion of these forces is defined by the connecting points of bodies i and j (see 
Fig. 4.7) where the actuator is installed. A vector l tJ connecting points P x 
and Pj is defined as 

lij = rfe + if b t -rje- if bj (4.7.3. 1) 

The actuator force f a acting on bodies i and j is given by 


fi l) = ±fJa ; /i j) = T/a/a 


(4. 7.3. 2) 
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where l a = The ± sign constitutes a pull and a push forces that are 

given by the actuators. The contribution of /i l) and / to the force vector 
F a on bodies i and j are 

F W> _ [/<«>, *f. y) x (4 7.3.3) 



Fig. 4.7 An Actuator Acting on Two Bodies 


4.7.4 Damping forces 

Dampers dissipate relative body motions by converting mechanical 
energy to dissipated forms such as heat. The damping force between bodies 
i and j at point Q t and Qj (Fig. 4.8) is found to be 


fd = d 


ip.j 


m 


ij**3 I 


where d is the damping coefficient and 


(4.7.4. 1) 




(4. 7. 4. 2) 


The damping forces acting on bodies i and j are 


/j ^ — -ffdld 


(4. 7. 4. 3) 
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in which Id = yj^j. The contribution of and to the force vector F<i 
on bodies t and j can be found from (4. 7.3.3). 



Fig. 4.8 A Damper Acting Between Two Bodies 
4.7.5 Spring Forces 

In mechanical systems, springs are often used to restore position 
equilibrium of two bodies. In Fig. 4.9, a spring is attached to two points S t 
and Sj of bodies i and j. The spring force is calculated by 


fs = k{ld-lo) (4.7.5. 1) 

where k is the spring coefficient, Id = |Sy — S t -| is the deformed length and 
lo is the undeformed length along the vector between two points S t and S ; . 
The spring forces acting on bodies i and j are 

/i l) = ±fjs ; / 3 (J) = Tfjs (4. 7.5.2) 

where l a — 4|4t. The contribution of /j 1 ^ and /j ^ to the force vector F a on 




66 


bodies i and j can be found from (4. 7. 3. 3). 



4.8 Concluding Remarks 

This chapter gives explicit mathematical expression for mechanical 
joints and forces pertaining to MBD systems. It is emphasized that the con- 
straint Jacobian matrix is obtained by extracting velocity vectors from the 
time differentiation of the constraint equations. Under such circumstances, 
each joint, which is represented by a different constraint Jacobian matrix, 
can be written separately without risk of confusion. From a programming 
standpoint, this development enables MBD software to possess modularity 
so that the equations of motion for multidisciplinary engineering problems 
can be automatically generated. 

The remaining issues regarding DAEs emphasize solution proce- 
dures. Chapter 5 reviews existing solution procedures, their advantages 
and disadvantages, and proposes two new constraint treatment schemes in 
conjunction with the two-stage staggered explicit-implicit algorithm to form 
a numerically robust solution procedure. 




CHAPTER V 


SOLUTION PROCEDURES FOR DAEs 

5.1 Introduction 

The equations of motion for MBD systems that are formulated in 
chapters 3 and 4 are characterized as DAEs. Since a closed form solution 
to DAEs cannot be found except for highly simplified problems, numerical 
methods must be applied in order to solve these highly nonlinear equations. 
Several existing numerical methods for solving DAEs are discussed in sec- 
tion 5.2. These numerical methods have primarily focused on the treatment 
of the constraint equations, which involved either constraint stabilization or 
constraint elimination. However, while these methods offer a varying degree 
of success, the lack of broadly applicable and robust numerical algorithms 
for solving DAEs remains as a challenging topic in the field of MBD systems. 
In this regard, two robust numerical methods that deal with both constraint 
stabilization and constraint elimination of DAEs are developed in sections 
5.3 and 5.4. In section 5.5, a numerical algorithm called two-stage stag- 
gered explicit-implicit procedure is developed to integrate translational and 
rotational motions separately. The stability criteria of this algorithm will 
be derived by linearizing the Euler equations. It is shown that the present 
algorithm not only prevents the instability but also maintains the explicit 
nature of the algorithm. 

5.2 Reviewing of Existing Solution Procedures 


In reviewing existing DAEs solution procedures, a numerical method 
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that is based on Gear’s backward difference algorithm has been first devel- 
oped to solve DAEs [20], With little success by using this algorithm, Gear 
and Petzold [21,22] have solved DAEs by differentiating the constraint equa- 
tions twice in time and augmenting these equations with the governing equa- 
tions of motion to form a combined set of second-order differential equations. 
If the augmented constraint equations are numerically integrated, however, 
constraint violations will generally occur because of accumulated integration 
errors. Several solution procedures that are based on the generalized coor- 
dinate partitioning technique [27,53], Baumgarte’s constraint stabilization 
technique [23,24], and the null space method [30-32,54] have been developed 
to overcome this key drawback. In the following sections we address these 
solution procedures in detail. 

5.2.1 Stiffiy-Stable Gear Method 

The earliest numerical algorithm that was used to solve DAEs is the 
so-called stiffly-stable Gear algorithm [20] that has been applied to some re- 
stricted differential-algebraic equations. This algorithm was used to form 
a set of augmented equations of motion by considering the algebraic con- 
straint equations as a special form of differential equations in which time 
derivatives of the variables do not appear. The equations of motion are then 
substituted into the backward difference formula and solved simultaneously 
with algebraic constraint equations that represent the kinematic joints of 
the system. The method starts with transforming DAEs into the following 
equations as 

Mv + B t A = F (5.2. 1.1) 


U = V 


(5.2. 1.2) 
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*( q,0 =0 

which can be expressed into the following matrix form 


(5.2. 1.3) 


f(y,y,0 = 


where 


M 

0 

o' 

(v I 


1 

td 

> 

0 

I 

0 

u 

+ 

“ v 

0 

0 

0_ 

UJ 

1 

* J 


G = 


g 


MOO 
0 10 
0 0 0 


F - B t A 


—v 

$ 


= Gy + g = 0 (5.2. 1.4) 


y = [q,v,A] T ,y = (q, v, A 


Solutions of (5.2. 1.4) may contain both high and low frequency components 
depending upon the driving term g in (5. 2. 1.4) and the eigenvalues of the 
system, the present system may become numerically stiff. Numerically stiff 
systems are characterized by having solutions dominated by low-frequency 
components. However, due to the presence of high-frequency components, 
the time step of the explicit numerical algorithms must be kept relatively 
small. This means that a large number of time steps is required to obtain 
accurate solutions. Consequently, schemes that damp out errors associated 
with high-frequency components are desirable. The Gear algorithms, due 
to their stiffly-stable characteristics for high-frequency ranges, are thus well 
suited for solving stiff DAEs with parabolic characteristics. 

The solution procedure starts with the Newton-Raphson algorithm, 
which is adopted to compute y, applied to (5. 2. 1.4) as 


f W Ay ( fc ) + f^ fc) Ay (/c) = -f (/c) 


(5.2. 1.5) 
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where k is the iteration cycle number. To obtain the relation between Ay^) 
and Ay^, the p th order Gear algorithm for (5. 2. 1.5) is employed and writ- 
ten as 

p-i 

y i+1 = 2Z(a jy *->) + 6 h f(y‘ +I ,y ,+1 ) (5.2.1.6) 

t=l 

where h is the time step, and ay and b are the coefficients that need to be 
determined depending upon the order of the algorithms one wishes to apply. 
For the k th and (/c + l) th iterations, (5. 2. 1.6) can be rewritten as 

p-i 

(y i+1 )‘ = (£(«yy i_, '))‘ + * f(y i+1 ,y’ +1 )‘ ( 5 . 2 . 1 . 7 ) 

*=1 

p-1 

(yf+i)fe+i _ (y^(ayy i_y ))* +1 + bh f(y v+1 ,y i+1 ) fe+1 (5.2. 1.8) 

t=l 

Since the summation term of (5.2. 1.7) and (5.2. 1.8) are only a function of 
the i th and previous time steps, they remain constant during the iteration. 
Subtracting (5. 2. 1.7) from (5. 2. 1.8) yields 

Ay (fc) = -!-Ay (fc) (5.2. 1.9) 

o a 

which gives the relation of Ay( fc ) and Ay^ fc \ Substituting (5. 2. 1.9) back 
into (5. 2. 1.5), we obtain 

(f» + Af 4 )“>Ay<*> = (g, + A G )l‘>Ay< fc > = -f‘*> (5.2.1.10) 

Upon solving (5.2.1.10), the update value of y and y can be obtained by 


y (*+i) = y (*) + Ay( fc ) 
y(c+i) = yW + -i-Ay(*> 


(5.2.1.11) 

(5.2.1.12) 


The drawbacks of this procedure are: first, it has expanded n -r rn DAEs into 
2n + m equations, thus for a large-scale system, it presents some inefficiency 
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solving y and y by using Newton-Raphson algorithm. Second, at starting 
time, there is no information available to determine the Lagrange multipliers 
A. An poor estimated A may cause the system to diverge. Third, due to the 
requirement of many time steps, the time discretization may have compound 
effects on the constraint equations, and eventually lead to useless drifted 
solutions. 


5.2.2 Direct Integration Method 

In the previous section, we have observed the difficulty of choosing 
Lagrange multipliers. To overcome this difficulty, Gear and Petzold (21,22 
purposed to convert the DAE system into a set of ordinary differential equa- 
tions by appending the second time derivative of the constraint equations 
to the state equations. Combining (3.4.19) and (3.4.20) with the governing 
equations of motion, the resulting constraint-augmented equation can be 
written in the following matrix form: 


M 

B 



(5.2.2. 1) 


where c = -(Bu + 2$ ut u + $«). The Lagrange multipliers A can be 
calculated by solving ii of the first row of (5.2.2. l) and using the resulting 
expression substituted into the second row of (5.2.2. l) so that 


BM _1 B r A = BM _1 F + c (5. 2. 2. 2) 

If the m x m matrix BM -1 B r is not singular, the acceleration vector ii 
can be computed by substituting the result of (5. 2. 2. 2) in to the first row of 
(5.2.2. l) to yield 

ii = M -1 [F - B r (BM - 1 B r ) - 1 (BM _1 F + c)] 


(5. 2. 2. 3) 
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Numerical algorithms are then employed to compute the generalized ve- 
locity and position vectors at the next time step. Note that, upon using 
this approaches, several difficulties may occur: First, during the process 
of simulation, BM -1 B r may become ill-conditioned, which degrades the 
numerical accuracy of (5. 2. 2. 2). Second, because numerical integration al- 
gorithms are used to compute ii, difficulty will occur for the reason that 
numerical time integration algorithms provide only an approximate solu- 
tion of the equations. During the time integration, the numerical errors 
may start to accumulate to the point that the constraint conditions are no 
longer satisfied to the desired accuracies. Since there is no numerical mech- 
anism to correct this defect, the solution may gradually diverge from the 
exact solution. This numerical difficulty is known as constraint violations. 
Third, the second time derivatives of the constraint equations do not nec- 
essarily represent the original algebraic constraint equations with fidelity in 
the case where nonlinear expressions are involved. 

To avoid aforementioned difficulties, several corrective methods have 
been proposed. These methods include the generalized coordinate partition- 
ing, Baumgarte constraint stabilization, the penalty constraint, and null 
space, which are reviewed in the sequel. 

5.2.3 Generalized Coordinate Partitioning Method 

The generalized coordinate partitioning method (GCPM) was first 
developed by Calahan [53]. Wehage and Haug [27] use this idea to reduce 
the system equations and determine independent coordinates from the con- 
straint equations. Their approaches are based on the fact that the n gener- 
alized coordinates of DAEs are not all independent. If the n coordinates can 
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be partitioned into m dependent and n — m independent coordinates, then 
the velocity and acceleration vectors can be partitioned into u d , u', u , and 
u l accordingly, where d, denotes the dependent generalized coordinates, and 
i denotes the independent generalized coordinates. The constraint equations 
(3.4.14) and its time derivatives (3.4.18) and (3.4.19) can be rewritten into 
the following forms: 

*(u d ,u < )=0 (5.2.3. 1) 

B'V = — B‘u* (5.2. 3. 2) 

Bii + Bti = [B d |B i ] | " t - j + Bu = 0 (5.2.3.3) 

Since B d has m full row rank, which indicates |B d | ^ 0, therefore the 
dependent acceleration vector is given by 

u d = -B d-1 (B‘u‘ +Bu) (5.2. 3. 4) 

The equations of motion (3.4.13) can be rewritten into the following parti- 
tioned form 



^ 

Premultiplying by B d on the first row of (5. 2. 3. 5) yields 

B d_T M d ti d + A = B^"^ 4 (5. 2. 3.6) 

Substituting A of (5. 2. 3. 6) into the second row of (5. 2. 3. 5) yields 

M‘u‘ + B ,T (B d_T F ti - B d ~ T M d u d ) = F‘ (5.2.3. 7) 

Replacement of u d in (5. 2. 3. 7) by (5. 2. 3. 4) leads to 

(M* + B^M^B*)!! 1 = F* - B' T F j - B* r M d B d_1 Bu (5. 2. 3. 8) 
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where B* = B d B*. Equation (5. 2. 3. 8) represents the reduced form of 
DAEs where a set of independent differential equations is given. Numeri- 
cally, these independent differential equations can be solved without violat- 
ing the constraint equations. From a computational point of view, (5. 2. 3. 8) 
is used to integrate the independent variables, whereas the dependent vari- 
ables are obtained by satisfying the constraint equations (5. 2. 3.1), (5. 2. 3. 2), 
and (5. 2. 3. 3) at each time step. During the process of solving independent 
acceleration vector, one may not under estimate the computational cost of 
factorizing the left hand side of the fully populated matrix (5. 2. 3. 8). A 
GCPM algorithm is stated as follows: 

(1) Given initial conditions u°, and u°. 

(2) Solve (5.2.2. 1) for u n , and A n . 

■* ,n 

(3) Specify (Compute) independent variable u l . 

*■ n-j- 1 * .rc-h 1 

(4) Integrate u* to obtain u 1 , and u l 

(5) Solve (5.2.3. 1) for u dU+1 by using Newton-Raphson method ; u 71 ^ 1 is 
obtained. 

(6) Solve (5. 2. 3. 2) for u d and u n+1 is found. 

(7) Go to step (2) until the required run time is reached. 

In step (5), a good estimate of is needed so that the Newton- 
Raphson iteration may converge within a few iterations. A reasonable ap- 
proximation of u d can be obtained by taking the Taylor series expansion up 
to the third terms as 

jn-f 1 J7i ,’j n h? *' , n , . 

u d = + hu d H u d (5. 2. 3. 9) 

Wehage and Haug [27] have developed an algorithm to identify inde- 
pendent and dependent generalized coordinates by using L-U factorization 
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to decompose the constraint Jacobian matrix B. This algorithm occasion- 
ally leads to poorly conditioned matrices. When this occurs, it is necessary 
to choose a new set of independent generalized coordinates by repeating the 
L-U factorization process. This strategy not only increases the computing 
time but also propagates integration errors. In recent years, many research 
groups have developed numerical techniques such as the singular values de- 
composition [28,29] (SVD) and the QR decomposition [55] to factorize the 
constraint Jacobian matrix. These techniques provide some marginal ad- 
vantages over L-U factorization, but the main idea remains basically the 
same. 

5.2.4 Baumgarte’s Constraint Violation Stabilization Method 

To stabilize the constraint violations that occur in solving (5.2.2. l), 
Baumgarte [23,24] proposed a constraint violation stabilization method. 
This method modifies the original constraint equations to form a set of re- 
laxation differential equations that has the capability to suppress the growth 
of error and achieve a stable response. In this method, Baumgarte replaces 
the second row of (5.2.2. 1) by the following constraint equations: 

$ + 2<*$ + /? 2 $ = 0 (5.2.4. 1) 

for holonomic constraints, and 

$ + 7$ = 0 (5. 2. 4. 2) 

for nonholonomic constraints where a, /? 2 and qr are arbitrary positive con- 
stants for numerical stability. Obviously, 2 a$, f3 2 &, and 7 $ are the terms 
used to stabilize the error committed by the violation of constraint equa- 
tions and their time derivatives. To study the behavior of the method, the 
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general solution of the first order ordinary differential equation for constant 
7 (5. 2. 4. 2) is expressed as 

(5. 2. 4. 3) 

where the constant, A;,, are determined from initial conditions. Note that 
is decaying as time t is progressing. For the second-order ordinary dif- 
ferential equation (5. 2. 4.1), the general solution for constant a, 0 is 

= k u e^ lt + fc 2 .e M2 V = 1, ...,m (5.2.4. 4) 

in which A; lt and k 2 { are integration constants that depend on the initial 
conditions, and 

fi i t 2 = — a ± \/ a 2 — /? 2 (5. 2. 4. 5) 

In order to make the solution to the constraint equation decay optimally, 
the critically damped (5. 2. 4. 4) requires that a = (3 and a > 0 so that 

= {ku + k 2t t)e~ at (5. 2.4. 6) 

Substituting (3.4.18) and (3.4.19) into (5. 2. 4.1) and replacing u from the 
equations of motion, the Lagrange multiplier for holonomic constraints 
yields the following expression 

BM -1 B r A = BM _1 F + c + 2aBu + (5. 2. 4. 7) 

When a and (3 are equal to zero, equation (5. 2. 4. 7) recovers the original 
second time derivatives of the constraint equations (5. 2. 4. 3) where the nu- 
merical solution may diverge from the exact solution. For nonzero a and 
/ 3 , the numerical solution oscillates about the exact solution. The mag- 
nitude and frequency of these oscillations depend on the values of a and 
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j3. Park and Chiou (1986, 1987) have shown that for some MBD problems 
Baumgarte’s technique suffers from the following drawbacks: 

(1) When BM -1 B t is ill-conditioned, the accuracy of Lagrange multipliers 
is considerably degraded. 

( 2 ) The errors committed in the constraint equations decay with one time 
constant regardless of the physical natural of dynamic problems. 

(3) Large values of a, /?, and 7 will cause the damping terms to dominate 
the numerical solution of the equations of motion, and make the system 
numerically stiff. 

Since a and (5 play a key role in this procedure, an analysis of 
this method is undertaken to obtain relationships between the coefficients 
a, f3, and the time stepsize h. Mathematically, one can approximate the 
nonlinear constraint equations by using Taylor’s series expansion. In such 
an expansion, the holonomic constraint equations at the (f + h) time step 
are truncated after the second derivative terms: 

k 2 „ 

$(* + h) = $(t) + h*{t) + —&{t) (5.2. 4.8) 

in which t is the current time, and h is the time stepsize used in the numerical 
integration process. Since the constraint conditions should be vanished at 
time step (t + h), (5. 2. 4. 8 ) becomes 

+ = 0 ( 5 . 2 . 4 . 9 ) 

h h z 

Comparing (5.2.4. 5) with (5.2.4. l), a and (3 can be expressed in terms of h 
as 

1 

(5.2.4.10) 
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Replacement of a and (3 in (5. 2. 4. 5) by (5.2.4.10) lead to 

M1,2 = r( _1± 0. * = (5.2.4.11) 

It is noted that (5.2.4.10) does not reach the critical damped as one has 
shown in (5. 2. 4. 6). But if a constant time step integration algorithm is 
used, this modified version of Baumgarte’s technique indeed damps out the 
constraint violations faster than an arbitrarily assigned constant value of a 
and j3. Recently, Bae and Yang [56] have developed a method to determine 
the optimal stabilization constants a and /? so that the constraint violations 
can be damped out efficiently. But the difficulties concerning the appearance 
of an ill-conditioned BM -1 B T remain unanswered. 

5.2.5 Null Space Method 

The null space method [30,54] and its variations [17,18,31,32] are 
alternative methods to deal with DAEs that adopt special numerical pro- 
cedures to eliminate system constraint forces. Hemami and YVeimer [54 
introduced a matrix method by considering the m dimensional subspace 
spanned by the rows of the constraint Jacobian matrix. Let C be the or- 
thogonal complement of B, and A T be a (n — m) x n matrix whose rows 
span the subspace C. By definition 

A T B T =0 ' (5.2.5. 1) 

Premultiplying (3.4.13) by A r and using the result given by (5.2.5. 1), the 
governing second-order differential equations become 

A T Mii = A t F (5.2.5. 2) 

Since the choice of A r is not unique, the reduced system equations are not 
uniquely determined. Recently, de Jalon et al. ] 17,181 have utilized this 
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concept to obtain the null space of the constraint Jacobian matrix by using 
natural coordinates. Their procedure starts with the assumption that n-m 
independent velocities u l can be selected as a projection of u which can be 
expressed as 

u‘‘ = Eu (5.2.5. 3) 


where E is the matrix which defines the linear combination. Combining 
(5. 2. 5. 3) with (3.4.18) yields 


B 

E 


u = 



thus if 


B 

E 


is not singular, then 


(5.2.5. 4) 


u = 



If = 0, we conclude that 


u = Au* 


(5. 2. 5. 5) 


(5.2.5. 6) 


where A is an n x (n — m) matrix whose column constitute a basis of the 
null space of B. Replacement of u in (3.4.18) by (5. 2. 5. 6) leads to 


Bu = BAu 1 = 0 


(5.2.5. 7) 


But u l yt 0, which implies 

BA = 0 


(5.2. 5. 8) 


Transposing of (5. 2. 5. 8) gives 


A r B r = 0 


(5. 2. 5.9) 
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Equation (5. 2. 5. 9) has achieved the goal of constructing the null space of the 
constraint Jacobian matrix so that constraint forces can be eliminated. Aug- 
menting (5. 2. 5. 2) and (3.4.19), the final system of equations of the present 
procedure can be rewritten in the following matrix form 

’ A r M 
B 

Note that (5.2.5.10) not only destroys the symmetry of the matrix but also 
violates the constraint conditions when time integration algorithms are used. 
To avoid these drawbacks, one can either adopt the Baumgarte constraint 
stabilization technique or reduce and express (5. 2. 5. 2) in terms of system 
independent variables by taking the time differential of (5. 2. 5. 6) as 

u = Au‘ + Au l (5.2.5.11) 

Substituting (5.2.5.11) into (5. 2. 5. 2) yields 

A r MAu‘ = A r F - A t MAu' (5.2.5.12) 

This expression eliminates the system constraint forces and achieves the goal 
of symmetrizing system equations without violating constraint equations. 
However, this approach suffers from two major drawbacks, the first one 
arising from the fact that if numerical integration algorithms are used to 
integrate the independent accelerations u\ the null space matrix A and its 
time derivative A need to be evaluated in advance so that the independent 
acceleration can be determined. Since the null space matrix A and its time 
derivative A are obtained by solving system dependent velocity and position, 
resolution of the constraint equations by using a costly Newton-Raphson 
iteration becomes unavoidable. The second drawback is that during the 


- -2 


(5.2.5.10) 
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process of time integration some positions can be reached which can cause 
the matrix in (5. 2. 5. 4) to become singular. The reason is that the chosen 
independent variables do not numerically represent uniquely the motion of 
all the possible mechanisms during the process of simulation. To be able to 
continue the integration process, a new matrix E must be chosen. 

Since these solution procedures suffer to some degree from com- 
putational difficulties, we have motivated to look for alternative solution 
procedures that overcome such difficulties. These new solution procedures 
which involve either constraint stabilization (penalty constraint stabilized 
technique) or constraint elimination (natural partitioning scheme) will be 
developed in the following sections to emphasize their numerical robustness 
and efficiency. 

5.3 Penalty Constraint Stabilization Technique 

In Baumgarte’s method, the objective is to minimize the error ac- 
cumulated in the constraint equations. In the penalty technique, instead 
of trying to eliminate the constraint violations, the errors being committed 
in constraint equations will be controlled. In other words, by applying the 
concept of proportional control of the constraint equations, the Lagrange 
multipliers are determined from the violation of the constraint equations as 
(Lanczos, 1949 [51]) 

$ 

A = — , 0 < c << 1 (5.3.1) 

e 

where e is the penalty coefficient. Substituting (5.3.1) into the equations of 
motion, the final equation becomes 

Mii + -B t $ = F 

e 


(5.3.2) 
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An important consideration in using penalty techniques is that we assume 
the constraint violations always exist. A reduction of the constraint vio- 
lations may be accomplished by decreasing the penalty coefficient which 
provides large number for corresponding Lagrange multipliers. As a result, 
the system equation (5.3.2) will be greatly stiffened which makes the in- 
tegration of the equations of motion become numerically difficult. On the 
other hand, if the penalty coefficient is increased, the errors propagate into 
the integration process and may lead to an unacceptable drifted solution. 
Furthermore, the penalty coefficient can not be a fixed constant. The rea- 
son is that during the process of integrating a constrained dynamic system, 
different constraint equations may require different penalty coefficients in 
order to stabilize different constraint violations. This motivates us to look 
for an alternative way to stabilize the constraint violations. 

To overcome the difficulties that have occurred in the previous tech- 
niques, a penalty constraint violation stabilization procedure [19,20] has 
been successfully introduced to correct the errors accumulated in the con- 
straint equations. This procedure is based on the observation that a time dif- 
ferential equation of penalty formula retains the characteristic of parabolic 
in time so that it is amenable to direct time integration and the constraint 
violations will decay according to the different time constant. To illustrate 
this procedure, the nonholonomic constraints will be treated first. The 
penalty technique for the nonholonomic constraints is expressed as 

eA = Bu + $(t) (5.3.3) 

Time differentiation of (5.3.3) yields 


eA = Bii -fi Bu -I- 


(5.3.4) 
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Substituting u from the equations of motion into (5.3.4) leads to 

eA + BM _ 1 B T A = Bu + BM _1 F + $ t (5.3.5) 

If A = ye ~' ,t , the eigenvalues of homogeneous part of (5.3.5) are 

( 7 fc I+ BM £ lB - T )y = 0 (5.3.6) 

where ( 7 fc , k = l,...,m) are the eigenvalues of (5.3.6). Since 7 *, dictates how 
the errors will decay with time, the different time constants will correct the 
errors committed in the constraint equations. This property overcomes the 
difficulty in Baumgarte’s method that errors that have been accumulated in 
the constraint equations can only decay with a fixed given time constants. 
For holonomic constraints, time derivative of (5.3.1) yields 

eA = Bu + (5.3.7) 


Integrating the equations of motion once to obtain the velocity vector, we 
get 


u n = h£ + 2 M~ 1 (F-B T A) n (5.3.8) 


where g is half of the time stepsize /i, and is the historical velocity vector 
that depends upon the applied numerical algorithms. Substituting \\ n into 
(5.3.7) yields 


eA" + ? BM- 1 B r A n = B n (^M -1 F + h£) n + $ t (5.3.9) 

Regardless of the nature of the constraints as shown in (5.3.5) and (5.3.9), 
the computation of Lagrange multipliers will not cause any numerical dif- 
ficulty even if BM _ 1 B T becomes ill-conditioned. Furthermore, this tech- 
nique provides two sets of differential equations for solving generalized coor- 
dinates and constraint forces. Since these two differential equations can be 



84 


solved separately, the software modularity of present procedure is enhanced. 
Several example problems have been examined and shown that the penalty 
constraint stabilization technique not only treat the constraints correctly 
but yields more robust solutions. 

In the present section, we have developed a procedure to minimize 
the constraint violations without introducing any artificial damping. In 
chapter 7, several example problems will be used to demonstrate the ro- 
bustness and efficiency of the present stabilized technique. An ultimately 
different approach to prevent constraint violations will be introduced in the 
next section by adopting the concept of the null space method. 

5.4 Natural Partitioning Scheme 

A partitioning scheme based on physical-coordinate variables is pre- 
sented in this section to systematically eliminate system constraint forces 
and yield the equations of motion of multibody dynamics systems in terms of 
their independent coordinates. Key features of the present scheme include: 
First, an explicit determination of the independent coordinates. Second, 
a parallel methodology to construct the null space matrix of the constraint 
Jacobian matrix is addressed if system topologies consist of a number of tree 
structures. For a system that contains closed-loops, a cut-joint technique 
is used so that the present scheme can be applied. Third, an easy incor- 
poration of the previously developed two-stage staggered explicit-implicit 
solution procedure. 

5.4.1 A Single Open Chain MBD System 

To demonstrate the present physical coordinates partitioning scheme 
for open loop systems, a three-dimensional triple-pendulum problem is cho- 
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sen (Fig. 5.1). 



mg 


Fig. 5.1 Example of Open Chain a MBD System: 
Three-Dimensional Triple Pendulum 

The constraint equations for this problem can be written as 

'[B n ] 0 0 

[B 21 ] [B 22 ] 0 < U2 > = 

0 [B3 2 j [B33] _ \ U3 J 

'[-/] [R a n] 0 0 0 0 1 ( ill ' 

[7] [R s21 ] [-1] (R 322 [ .0 0 ln 2 > 

.0 0 [I] [R332] [-1] [R^as] J l u 3 , 


= 0 (5.4. 1.1) 
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where the pendulum bodies are connected by three spherical joints and R s 
are the function of rotational operators and position vectors from the center 
of mass of each body to the position of their connecting joints. To obtain 
the necessary projection matrix A, we start with the first row of (5. 4. 1.1): 


BnUi = [ -I , R aU ] Ui = 0 (5.4. 1.2) 


that can be partitioned into 

[Bf 1 |B\ 1 ||"j}=0 (5.4. 1.3) 

or 

B d 1 u d + B\ l u\ = 0 (5.4. 1.4) 

where Bf x = —I, B l n = R s n, d represents the dependent coordinates and 
i represents the independent coordinates. Since |Bj X | ^ 0, the dependent 
velocity components of first body can be calculated as 


Uj — Bjj BjjUj — PiUj 


(5. 4. 1.5) 


where Pj = — Bjj B' n = R a n- The velocity vector of first body lit can 
be written in terms of independent velocities li ^ as 


Ui = 


u l 


P j ) ui = qx 


(5.4. 1.6) 


where Q x 


Pi 


. Likewise, B 22 of the second row of (5.4. 1.1) can be 


partitioned into 


B 2 iUi + 822^2 + 622^2 = d 


(5.4.1. 7) 


or 


U 2 — — B 2 2 (^ 21^1 + ^22^2) 


(5.4. 1.8) 
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for |B 22 | 7^ 0 - Substituting ( 5 . 4 . 1.6) into ( 5 . 4 . 1.8) yields 


*2 = -B22 1 (B2iQ 1 u‘i + B‘ 22 u l 2 ) = Riu\ + R 2 u 2 


(5.4. 1.9) 


where Ri = — B 22 B21Q1 = B21QH and R2 = — B 22 B 22 = B 22 . The 

velocity vector of the second body, U2, can be expressed in terms of the 
independent velocities, and u 2 , els 


u 2 


• d. 

^ } U ? I = 

U, 


Ri R2 
0 I 


u i 

111 


[S1IS5 


U 1 

Ur 


(5.4.1.10) 


where Si = 


^ ^ and S2 = ( ) . Applying the same procedure to the 


third row of ( 5 . 4 . l.l), u 3 can be expressed as 

U3 = — B 3 3 (B 32 U 2 + B33U3) = — B3 3 [B32(SiUi + S 2 U 2 ) + B33U3] 

= v^i + V2U2 +V3U3 (5.4.1.11) 


where Vi = — B33 B32S1 = B32S1, V2 = — B33 B32S2 — B32S2, V 2 


B33 ^33 = B33, and U3 can be written in terms of u\, u 2 , and 113 as 


u 3 



'Vi 

V 2 

V 3 ‘ 

Uj _ 

0 

0 

I 



where W x 


^o)’ W2 " \° 2 J ’ ““ 

( 5 . 4 . 1 . 6 ), ( 5 . 4 . 1 . 10 ), and ( 5 . 4 . 1 . 12 ), we construct the physical velocities u in 


I U ‘ 

W 1 1 W->| W 3 j l uf, 

1 4 

(5.4.1.12) 

Vo \ 

and W 3 =1 _* I . Combining 


terms of u' as 



Qi 

Si 


0 

S2 


0 

0 


or 


W, Wo W 3 
u — Au* 



(5.4.1.13) 


(5.4.1.14) 
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where A is the null space of the constraint Jacobian matrix that has been 
exploited in the previous section. Note that in the process of forming A, the 
inversion of the dependent matrices can be obtained analytically as opposed 
to the generalized coordinate partitioning scheme in which the inversion of 
the dependent matrices has to be carried out numerically. The scheme for 
constructing A provides a guideline to deal with MBD systems containing 
different topologies such as multiple open kinematic links and closed kine- 
matic loops as discussed in the sequel. 

5.4.2 A Multiple Open Chain System 

If the MBD systems have more than one branch as shown in Fig. 5.2, 
the present scheme lends itself to multiprocessor computers. This property 
can be demonstrated by the following simple MBD system for which the 
constraint equations are given by 


riBnl 

0 

0 

0 

0 


• 

Ul 

[Bai] 

[B 22 ] 

0 

0 

0 


u 2 

0 

[b 32 ] 

[B 33 ] 

0 

0 

< 

U 3 

|B«] 

0 

0 

|b 44 ] 

0 

i 

U4 

0 

0 

0 

[B54] 

(B 55 ] . 


, U5 ) 


Applying the scheme used in section 5.4.1 to both chains, the null space of 
the constraint Jacobian matrix A is decided as 


f ^ 


‘Qi 

0 

0 

0 

0 ' 


{ • t \ 

U 1 

U 2 


Si 

S 2 

0 

0 

0 


U2 

u 3 

y = 

w x 

W 2 

W 3 

0 

0 

i 

U3 > 

u 4 


Yx 

0 

0 

Y 4 

0 


^4 

, u 5 , 


. Zi 

0 

0 

z 4 

Z 5 . 


U l 5 J 


(5.4. 2. 2) 
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Fig. 5.2 Example of MBD Systems with Multiple Branches 

Note that, in the physical coordinates partitioning scheme, once the 
first row of (5.4.2.2) is constructed, the second and fourth row of (5. 4. 2. 2) 
can be constructed simultaneously according to the given Q i . Again, if the 
first, second, and fourth rows of (5. 4. 2. 2) are found, the third and fifth rows 
of (5. 4. 2. 2) can be obtained according to their dependent branches respectly. 
Since MBD systems are the systems that include many kinematic loops, it is 
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natural to utilize this development in a multiprocessor computer to compute 
the null space (at each branch) of the constraint Jacobian matrix. 

5.4.3 A Closed-Loop MBD System 

When the systems have one or more closed loops, difficulties arise 
in constructing the null space of the constraint Jacobian matrix as one can 
see by examining the following three body crank-slider problem (Fig. 5.3). 



Fig. 5.3 Example of a Closed-Loop MBD System: 
The Crank-Slider Mechanism 


The constraint equations for this problem are given by 


'[Bn] ' 

0 

0 


/ • \ 

(B 2 1 ] 

[B 22 ] 

0 



0 

[b 32 ] 

[b 33 ] 


1 

0 

0 

!B43] 


l u 3 J 


= 0 


( 5 . 4 . 3 . 1 ) 


It is obvious that joint 1 and 4 conflict in determining the null space of 


(5.4.3. 1) according to the preceding scheme. Fortunely, there is an elegant 



91 


way of overcoming that difficulty. The technique relys on “cut joints”, which 
means cut the joints that are necessary to force the system topology to 
become open loop so that the existing solution procedure could be adopted. 
This technique is accomplished by partitioning (5.4.3. 1) into the following 
form 


' [Bn] 

P21] 

0 

. 0 0 |B„1 J 

or 

B o ii = 0, B c ii — 0 (5. 4. 3. 3) 

where B 0 represents the open loop constraint Jacobian matrix, and B c rep- 
resents the remaining constraint Jacobian matrix after the joints have been 
cut. Performing the physical coordinates partitioning scheme to construct 
the null space of B a as 

B 0 A 0 = 0 ; = 0 (5. 4. 3. 4) 

Performing the algebraic calculations as in section 3.1 yields the equations 
of motion for a closed-loop MBD system as 

Mu + B^A 0 + BjA c = F (5.4.3. 5) 

Premultiplying A^ to the above equation yields 

AjMu + = AjF (S.4.3.6) 


0 

[B 22 ] 

B3 2 j 


0 

0 

[B33] 



!»:}* 


(5. 4. 3. 2) 


or 


M n ii + B^A C = F 


( 5 . 4 . 3 . 7 ) 
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that are subjected to the constraints 

B c ii = 0 (5. 4. 3. 8) 

where 

M n = AjM 

Bl = AjBj (5. 4. 3. 9) 

F„ = AjF 

Equations (5. 4. 3. 7) and (5. 4. 3. 8) can be solved either by employing the 
penalty constraint violation procedure that has previously developed or by 
constructing the null space for the new equations of motion. 

5.5 Explicit-Implicit Solution Procedures 

In sections 5.3 and 5.4, two schemes have been developed to treat 
constraints efficiently. The remaining task for the numerical solution of the 
equations of motion of MBD systems is the computation of the general- 
ized coordinates. A solution procedure called two-stage staggered explicit- 
implicit procedure [43] has been developed to integrate the governing equa- 
tions of motion. This procedure based on the partitioned solution procedure 
hats been used to partition the governing equations of motion into transla- 
tional and rotational components. Two numerical algorithms are used to 
integrate the generalized coordinates and constraint forces of the penalty 
constraint stabilization technique, and the generalized coordinates and in- 
dependent coordinates of the natural partitioning scheme. The following 
sections describes this procedure in more detail. 

5.5.1 Partitioning the Governing Equations of Motion 

By using the penalty constraint stabilization technique, the discrete 
equations of motion for MBD system derived in section 3.4 may be expressed 
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as 


Mu + B 1 A = F 


( 5 . 5 . 1 . 1 ) 

eA = Bu + ®t 

It is noted that the integration of angular velocity does not lead to angular 
orientations, which are obtained by integrating a separate set of kinematic 
equations involving angular velocities. This motivates us to partition the 
generalized coordinates u into translational velocity vector, r and angular 
velocity vector oj into the following forms 


*■(«)■ *■(*) 


( 5 . 5 . 1 . 2 ) 


so that the desired angular orientations can be obtained by solving a set 
of kinematic equations. Since the translation and the rotation components 
are decoupled in the mass matrix, the equations of motion can be further 
partitioned into 

>r 


Mr 0 

0 M, 


uj 


+ 


T ? A 


F r 

F,„ 


( 5 . 5 . 1 . 3 ) 


where the subscripts (r, cj) refer to translational and rotational motion re- 
spectively. The translational and rotational parts of the constraint .Jacobian 
matrix are given by B r and B w . Note that (r,u>) can be partitioned into 
body-by-body degrees of freedom as 


r = [r 1 ,r 2 , ,r n6 ] r 


oj = [a > l ,oj 2 , , 


, OJ 


nb i T 


( 5 . 5 . 1. -1 1 


where nb is the total number of bodies in the system and r ! and oj 1 are the 
translational and rotational velocity vectors for the i-th body where 


f* _ [ft f« ft \T 
r — l r l i ~2’ r 3l 


oj 1 = [u;‘ 1 ,u;^,u;^] r 


( 5 . 5 . 1 . 5 ) 
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5.5.2 Two-Stage Staggered Explicit-Implicit Procedure 

Since DAEs are not differential equations, a special study is needed 
to select robust and stable algorithms that yields accurate solutions. In this 
regard, studies must be conducted in order to analyze different numerical 
algorithms with their stability and accuracy characteristics. In structural 
dynamics, the most widely used numerical formula of solving the discrete 
equations of motion is the central difference formula, the half-station version 
of which may be written: 


u 


u 


n + k 


= u 


+ hii T 


n+1 = u n + hi i re+ = 


(5.5.2. 1) 


where h is the time stepsize. This numerical algorithm is an explicit inte- 
gration algorithm with second order accuracy. If u max is the highest in- 
stantaneous frequency, the stability condition of central difference formula 
is uj max h < 2 which imposes the time stepsize limitation. This algorithm is 
attractive to parallel computations because of its robustness and explicit na- 
ture, therefore the application of present algorithm to MBD systems needs 
to be evaluated carefully. A dynamic torque-free top problem will be demon- 
strated here to investigate the drawback of the central difference algorithm 
if the equations of motion are a function of velocity. The torque-free top 
problem is governed by the equation 


Juj + cu x Jiu — 0 (5. 5. 2. 2) 

where J is the inertia tensor. If the central difference formula is used to 
integrate this equation of motion, the angular velocity u is obtained at the 
half time step via the integration formula, while the angular acceleration ^ 
are obtained at the full time step via the governing equation. Since w at 
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the full time step are not available, the angular acceleration Co at the full 
time step cannot be found. In order to continue the integration process, two 
approximations have been studied: 

(1) to n ~ U} n ~5 

(2) 0J n « |(w n ~5 + tu n+ ^) 

In order to assess the stability of these approximations, it is nec- 
essary to linearize the equtions of motion so that the stability criteria can 
be established. To linearize of (5. 5. 2. 2), we recall the rotational operator 
(2.3.1) and rewrite it into the following matrix form 

R(n, <f>) = n T n + cos <£( I - n r n) - sin <f>h (5. 5. 2. 3) 


The series for trigonometric functions sine and cosine are 


<f> 3 4> 5 

sm <£ = <£- — + — - + 


(5. 5.2. 4) 


e <t > 4 

cos <£=!- — + — - + 


(5. 5.2.5) 

Replacement of sin^> and cos <f> in (5. 5. 2. 3) by (5. 5. 2. 4) and (5. 5. 2. 5) lead to 


4> 


<j > 3 -T < t >4 


R(n,4>) =I 3x3 +— (n T n -I)-— n (n n - /) 


3! 


4! 


(5.5. 2.6) 


Next we observe the relations of the skew-symmetric matrix h that 


and 


(n T f = -h T , (ii T ) 4 = -(A T ) 2 


(n T ) 5 = n r , (h T ) 6 = (d T ) 2 , 


n T n — I = (n r ) 2 = — (n T ) 4 = (n T ) 6 = • • 


(5.5.2.T) 

(5.5. 2.8) 

(5. 5. 2.9) 



96 


Making the use of (5. 5. 2. 7-9), (5. 5. 2. 6) becomes 

R(n ,4>) = I 3x3 + <j>h T + ^-(n T ) 2 + |j-(» T ) 3 + ••• (5.5.2.10) 

or 

R(n ,<f>) = = e** (5.5.2.11) 

where 0 — <f>h . Since we have linearized the coordinate transformation 
matrix, the relation between b n and b n ~ 1 " 1 with an rotational angle 0 can 
be written into 

b n+1 = /V = e^R n e = R n+1 e (5.5.2.12) 

or 

R n+1 = R n (5.5.2.13) 

To establish the relationships between the angular velocity and the linearized 
parameters 0 , the angular velocity and the coordinate transformation matrix 
are related by 

<j n+1 = -R n+1 R n+i:r = — -j-(e^ R n )(R nT e^) (5.5.2.14) 

at 

Carrying out the time derivation of (5.5.2.14) leads to 

U n+X = 6 + J ) bJ n e 9 (5.5.2.15) 

n T 

Approximating e u of (5.5.2.15) by the first two terms of (5.5.2.10) which 
is given by 

a T 

e —1 — 0 (5.5.2.16) 

and substituting (5.5.2.16) back into (5.5.2.15) lead to 

<y * +1 « e + u n + Cj n e + e T u n 


(5.5.2.17) 
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or 

u n+1 & 0 + u n + CM (5.5.2.18) 

which can be expressed in vector form as 

w n+1 =u ) n + 6 + u n x9 (5.5.2.19) 

Taking the time differentiation of (5.5.2.19), the linearization of the angular 
acceleration u> n+1 is obtained as 


u> n+1 = u n + e + u n x 9 + u n x 0 (5.5.2.20) 

Replacement oftI; n+1 and in (5. 5. 2. 2) by (5.5.2.20) and (5.5.2.19) lead 

to 


Jw n+1 + u n+l x J • w n+1 =Ju n +u n xJ.u n + 

JO + [Jw n - (J^) + u n J\0 + [Jw n - (J^)w n + u n Ju n }0 = 0 (5.5.2.21) 

by neglecting any of the two linearized parameters product terms. The final 
linearized equations of motion can be expressed as 

JO + BO + K0 = 0 (5.5.2.22) 

where the gyroscopic damping D and the centrifugal force K are given by 
the following matrix forms: 

D = Ju n - (jZ/ 1 ) + u n J (5.5.2.23) 

K = Jw 1 - (J^)*» + w n Jw n (5.5.2.21) 

The central difference formula in (5.5.2. 1) can be algebraically transformed 
to 

0 n+l - 20 n + 0 n ~ l 


(5.5.2.25) 



■ 
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with the first velocity approximation 

•n • n — 5 

0^0 2 

in which 

n _x e n -e n ~ 1 

0 ~ — h 

Substituting (5.5.2.25-27) into (5.5.2.22) yields 

J(0" +1 - 20 n + 0 n ~ l ) + hB{0 n - 0 n ~ x ) + h 2 K0 n = 0 (5.5.2.28) 

The computational stability of this approach can be assessed by seeking a 
nontrivial solution of 

0 n + x = s 0 n = s 2 0 n ~ 1 (5.5.2.29) 

such that 

\s\ < 1 (5.5.2.30) 

for stability. Substituting (5.5.2.29) into (5.5.2.28) yields the following char- 
acteristic equation 


(5.5.2.26) 


(5.5.2.27) 


|J(s - l) 2 + hB{s - 1) + h 2 Ks\ = 0 (5.5.2.31) 

In order to examine the stability requirement on the characteristic equation, 
one transforms |s| < 1 onto the entire left-hand plane of the c-domain via 

5 = Ltl (5.5.2.32) 

1 — 2 

so that 

|J( T ^-) 2 + /,D( T ^) + /, 2 K(i^)| = 0 


(5.5.2.33) 
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It is noted that the stability criteria is often decided by the gyroscopic 
damping term rather that the centrifugal force term. To simplify the stabil- 
ity analysis, we set Ji = — J 3 = 1 and replace them into (5.5.2.23) and 

(5.5.2.24) to obtain 

D = £> (5.5.2.34) 

K =<i (5.5.2.35) 

Further simplification can be made if we substitute J\ = J 2 = Jz = 1 
back into the original governing equation (5. 5. 2. 2) to obtain u — 0 which 
implies K = 0. Substituting this expression and (5.5.2.34) into (5.5.2.33), 
the determination of (5.5.2.32) becomes 

\4z 2 I + 2hz(l — z)£j\ = 0 (5.5.2.36) 

Expanding (5.5.2.36), the resulting polynominal equation expressed in terms 
of z as 

[4-t-/i 2 (w 2 +w 2 +w|)]z 2 — 2/i 2 (tjj -fw 2 'h u; 3 ) 2 +^ 2 ( u; i +W 2 +W 3 ) =0 (5.5.2.37) 

The Routh-Hurwitz criterion asserts that for stability all the coefficients 
in (5.5.2.37) must be positive. Since the coefficient of the second term of 
(5.5.2.37) is negative, we conclude that the present velocity approximation 
makes the central difference algorithm numerically unstable. 

The second velocity approximation is 

6 n = -(0 n+ 5 + 0 n ~2) (5.5.2.38) 

or 

e n = \{0 n+l -r- 1 ) (5.5.2.39) 

2 h 
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By applying the same stability procedure, the z-polynominal equation is 
given by 

4z 2 + h 2 (u> 2 + Wj + wf) = 0 (5.5.2.40) 


Again, by using Routh-Hurwitz criterion, this algorithm is proven to be 
unconditionally stable. Thus we conclude that once the present algorithm 
is employed, a stable and accurate solution is obtained. However, since 
w n+ ^ is not available as part of integration process, the robustness and 
explicit nature of the algorithm have been lost. This has motivated us to 
develop the following two-stage staggered explicit algorithm which prevents 
the instability of the first velocity approximation while circumventing the 
unavailability of the second velocity approximation. 

In the two-stage staggered algorithm, the computational sequences 
have been divided into the following steps if 0 , 0 , and 0 are known: 


(1) d n+k - = e n * +hb n 

( 2 ) 0 n+ * = 0 n ~~* + h0 n 

(3) 9 + 2 is obtained by using (5.5.2.21) 

( 4 ) 0 n+1 = 9 U + h'9 n+h 

(5) 0 n+1 = 0 n + h0 n+h ‘ 

(6) <> n+1 is obtained by using (5.5.2.21) 

In this regard, we compute 0 by 

• rc • n— 1 . •:«- k 

9=9 + h9 - 

0 n+ 5 = 0 n ~5 + h0 n 


( 5 . 5 . 2 . 41 ) 


Following a similar step the same as in the previous cases, the characteristic 


equation for the two-stage algorithm becomes 


(s 4 - 2s 2 + 1)1 + hu{s 3 - s)\ = 0 


(5.5.2.42) 
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Expanding and transforming (5.5.2.42), the resulting z-polynomial equation 
is given by 


h 2 Auz 4 + 2(8 — h 2 Au)z 2 + h 2 Au — 0 (5.5.2.43) 


where Au> = + wf. Equation (5.5.2.43) implies that the two-stage 

algorithm is stable if the stepsize remains in the range of 


h < 


2y/2 

(c^i + wf + w|)£ 


(5.5.2.44) 


The foregoing linearized stability analysis for the torque-free top problem 
confirms that the two-stage staggered algorithm not only provides a stable 
and accurate solution but maintains the explicit nature of the algorithm. 

For MBD systems, several example problems have been studied. 
A procedure called two-stage explicit-implicit procedure, which implements 
the previous development, has been used to integrate the governing equa- 
tions of motion. The two-stage explicit-implicit procedure uses the explicit 
central difference algorithm to integrate the translational and rotational 
velocities and the translational displacements. An implicit numerical algo- 
rithm is used to integrate the Euler parameters by imposing the kinematic 
relation between the angular velocity and the Euler parameters and their 
time derivatives. The following sections describe these procedures in more 
detail. 


5.5.3 Update of Translational and Angular Velocities 

Since (5. 5. 1.3) provides two sets of uncoupled differential equations, 
the translational and rotational acceleration vectors at n-step can be explic- 
itly computed assuming r n , u n , and A n are known. The translational and 
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angular velocity are integrated by the central-difference formula, namely 

T n+ 2 =i"-k + hr n 

U} n+ 2 = (jj n ~2 + h(l} n 

where the translational displacement r" + 5 is updated by using 


(5.5.3. 1) 


r n+§ _ r n ~ 2 + hr n 


(5. 5. 3. 2) 


Note that the updating of the next half step results by simply changing the 
index n + | ton + l,nton+^, and n - \ to n as the integration proceed. 

5.5.4 Update of Euler Parameters 

As indicated in section 2.8, the Euler parameters can be obtained by 
making the use of the kinematic relation (2.8.7). In the present algorithm, 
the Euler parameters are integrated by the following equation 

q n+1 = q n + hq n+k - (5.5.4. 1) 

Substituting (2.8.7) into the above equations yields 

q 1,1+1 = q n + fcA(w n+ ±)q n+ * (5. 5. 4. 2) 

The updated q n+1 needs to satisfy the constraint equation (2.6.3). Among 
several possibilities the approximation 

q n+ s = -(q n + q n+1 ) (5.5.1.3) 

2 

has been found to give the most accurate result. Replacement of q n+ ? in 
(5. 5.4. 2) by (5. 5.4.3) yields 

(I - ^(u.“ + i))q” +1 = (I + j A(u," + h)q" 


(5. 5. 4. 4) 
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The skew-symmetric matrices on the left hand side of (5. 5. 4. 4) can be ex- 
plicitly inverted via the formula 


'1 

—a 

-b 

— c 

-l 

1 

a 

b 

c 

a 

1 

c 

—b 

1 

—a 

1 

— c 

b 

b 

— c 

1 

a 

1 + a 2 + b 2 + c 2 

-b 

c 

1 

— a 

_ c 

b 

—a 

1 _ 


— c 

- -b 

a 

1 


(5. 5.4.5) 

Hence, the solution of q n+1 can be obtained without solving a set of linear 
equations as 

q " +1 = il I+ 5 A ( w " +i )l[I + ^(-“ + b]q" (5. 5.4.6) 

in which A = ~ (1 + uf + + w|). Once q” +1 is known, the coordinate 

transformation matrix R may be updated via (2.6.6) which relates the body- 
fixed basis of each body to the inertia basis. Note that for the two-stage 
staggered algorithm, the Euler parameters at n + i can be computed by 
shifting the index n + 1 to n + | and n + | to n as 

q “ + ' = il I+ (S.5.4.7) 


5.5.5 Update of Constraint Forces 

To compute the constraint forces for the hoionomic constraints 
(5.3.9), one integrates u n ^~^ and A n+ - with the following mid-point im- 
plicit numerical algorithm: 


=U n +QU n + k - 


A n+ 2 = A n + ^A 


n+ ; 


(5.5.5. 1) 


where q = |. When ii n+ - is substituted by the equations of motion, the 
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constraint equations for holonomic constraints becomes 
(el + 0 2 BM -1 B T ) n+ 3 A n+ ^ = eA n + gB n+ ^u n + 

(5. 5. 5. 2) 

The present two-stage explicit-implicit staggered procedure given by (5.5.3. 1- 
2), (5. 5. 4. 7), and (5. 5. 5. 2) constitutes a complete solution procedure for 
MBD systems that undergo large translations and rotations. The proce- 
dure of solving each quantity is given in the next section. 

5.5.6 Penalty Constraint Stabilization Technique Implementation 

[1] Initialize r,r,u;, and q. 

[2] Compute F at n step. 

[3] Compute initial A 

[4] Compute u = [r,u>] at n step. 

[5] Update translational and angular velocity r,u at n + | step by using 
(5.5.3. 1). 

[6] Compute q n+ 5 with (5. 5. 4. 7). 

[7] Compute A at n + j step. 

[8] Repeat [2] and [4] at n + | step. 

[9] Update r n+1 with (5. 5. 3. 2). 

[10] Repeat [5] and [6] at n + 1 step. 

[ 1 1] Extrapolate A to n + 1. 

[12] Go to [2] for next time step. 

5.5.7 Natural Partitioning Scheme Implementation 

[1] Initialize r,r,u>, and q. 

[2] Integrate r n+ 2 by r n+ 2 = r n- 2 + hr n 
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[3] Form A at n step. 

[4] Compute A r F — A T MAu t at n step. 

[5] Solve for n l at n step. 

[6] Compute u at n step. 

[7] Integrate u by using 



[8] Compute q n+ ^ according to (5. 3. 4. 7). 

[9] Repeat [2] at n + 1 step. 

[10] Repeat [3]-[6] at n + ~ step. 

[11] Repeat [7] and [8] at n + 1 to complete the step. 

[12] Go to [2] for next time step. 

5.6 Concluding Remarks 

The existing solution procedures of DAEs can be characterized into 
three categories: (l) treatment of DAEs as a special type of differential 
equation; (2) stabilization of the constraint violations by adopting special 
techniques; (3) construction of the null space of the constraint Jacobian 
matrix by employing appropriate matrix algorithms. So far, the relative 
performance of these methods have been measured largely in terms of a 
sequential computational context. However, new parallel computers are ex- 
pected to have significant influence on the algorithm development for the 
large-scale MBD systems and real-time simulation. To address this issue, 
two schemes, viz., constraint stabilization and constraint elimination, have 
been developed with parallel computation in mind. In the constraint stabi- 
lization scheme, the constraint violations that occur during the time integra- 
tion process are stabilized by adopting the rate form of the penalty scheme. 
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On the other hand, the constraint elimination scheme uses a decomposi- 
tion of the physical coordinates that manifests dependent and independent 
variables so that the null space of the constraint Jacobian matrix can be 
constructed explicitly. Because both schemes convert DAEs into special 
types of differential equations, an algorithm based on the explicit central 
difference formula and implicit mid-point rule was adopted to integrate the 
• equations of motion and their constraint forces or independent variables. 
Several example problems are used later to demonstrate the robustness and 
efficiency of this algorithm. 

After DAEs are solved, it is natural to search for a method that 
can dramatically reduce computer run-time. Again with parallel comput- 
ers in mind, this goal can be achieved by adopting a Schur complement 
based parallel preconditioned conjugate gradient numerical algorithm that 
determines: (l) system acceleration vectors and constraint forces for the 
constraint stabilization scheme; (2) system acceleration vectors and inde- 
pendent acceleration vectors for the constraint elimination technique. These 
details are covered in the next chapter. 



CHAPTER VI 


PARALLEL IMPLEMENTATION OF THE 
GOVERNING EQUATIONS OF MOTION 

6.1 Introduction 

During the past two decades, multiprocessing computer architec- 
tures have undergone progressive development because of the increasing 
availability of low cost multiprocessors. New parallel computers consist- 
ing of tens, hundreds, and even thousands of processors, have motivated 
the design of parallel algorithms and promised to have a profound impact 
on numerical simulation capabilities in many engineering disciplines. Some 
computer programs that run well on conventional sequential computers do 
not necessarily transformed to programs that efficiently harness the capa- 
bilities of parallel computers. This is particular true for massively parallel 
architecture. Conversely, algorithms that are less efficient in sequential com- 
puters may reveal an inherent parallelism that makes them attractive for 
parallel computers. 

Since an MBD system may consist of hundreds or even thousands 
of linked bodies, numerical solutions of such highly nonlinear systems may 
consume a large amount of CPU time. The ultimate purpose of real-time 
simulation has motivated the development of efficient parallel algorithms 
on existing parallel computers. The issues that directly pertain to parallel 
computations of MBD systems include: generation of the system equations 
of motion, incorporation or elimination of constraint forces, integration of 
generalized coordinates, and interpretation of the simulation results. As sug- 
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gested in previous chapters, the governing DAEs can be generated concur- 
rently without any particular difficulty. However, the solution process that 
involves constraint stabilization or constraint elimination may introduce dif- 
ficulties in the parallel determination of system generalized coordinates. In 
the present chapter a parallel solution method is proposed that computes 
directly the system constraint forces and generalized accelerations or system 
independent and generalized accelerations at the elementary body-by-body 
level. Once these quantities are known, the physical coordinates of each 
body can be computed independently and in parallel. 

This chapter is organized as follows: sections 6.2 and 6.3 review 
two MBD equations that were derived in previous chapters. By rewrit- 
ing these equations into body-by-body level, the governing equations can 
be transformed into a so-called Arrowhead matrix. The advantage of this 
matrix structure is that the generalized coordinates of system can be inte- 
grated simultaneously by using a previously developed two-stage staggered 
explicit-implicit algorithm. Section 6.4 discusses a parallel solution algo- 
rithm based on the conjugate gradient (CG) numerical algorithm, which is 
then applied to these arrowhead system equations. Finally, in section 6.5 
two CG preconditioners are introduced to improve the convergent rate ol 
the conjugate gradient numerical algorithm. 

6.2 Parallel Implementation of Penalty Constraint 

Stabilization Technique 

In order to obtain the optimum performance of the solution pro- 
cedure, a procedure is developed to utilize parallel processors for solving 
constraint forces and improving the computations involving BM B . In 
the two-stage explicit-implicit procedure, the numerical solution of MBD 
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systems is obtained by combining two modules: a generalized coordinates 
integrator and a Lagrange multiplier solver. The Lagrange multiplier solver 
integrates the constraint equations with the mid-point rule, namely 


= u n + 0 u n+ * 

A n+ 2 = A B + gA" + * 


( 6 . 2 . 1 ) 


where q is equal to the half of the time step k. When ii n+? is substituted 
from the equations of motion, the Lagrange multiplier solver is obtained. 
From the parallel implementation viewpoint, this solution scheme is not 
attractive on account of its complexity in computing BM _1 B r . An alter- 
native scheme is developed to compute the Lagrange multipliers more effi- 
ciently. The scheme uses previous derivations without substituting ii n+5 , 
which can be replaced by the governing equations of motion, into the La- 
grange multiplier solver. This leads to the following equations 

Mu n+ 5 + (B r A) n+ 2 = F n+ * 

(Bu) n+ * - 4A n+ * = -4A" - -B" + ^u n 

e 2 s 2 8 

or in matrix form 



(6.2.3) 


where E = -^1, c = F n+ *, d = -4r A n - jB n+ 5u". We partition M, u, 
and B as 



0 

0 

... 

£(1,0 

... 

£(l ,nj) " 


Ml 


'* Cl 

0 

( 2,2) 

0 

... 

£(2,1) 

... 

B (2,nj) 


&•} 


C‘2 

0 

0 

... 

Af(n6,n6) 

£(nt,l) 


B (nb,nj) 

< 

U n 

y = ^ 

c n ■ 

£(1.1) 

B ( 1,2) 

... 

B (\,nb) 

Ei 

0 

0 


Ai 


d i 

... 

... 

... 

... 

0 

. . . 

0 





- B {nj, 1) 

B (nj, 2) 

... 

B {nj,nb) 

0 

0 

E nJ . 


< ^nj j 

i 

. d -2 , 


V a 2 > 

( 6 - 2 . 1 ) 



no 


where nb and nj denote the total number of bodies and joints in the sys- 
tem respectively. Note the arrowhead matrix structure that appears in the 
acceleration equations can be written as 


nj 

Mj Uj + ^ = i 3 1 , 

i=1 (6.2.5) 

nb 

)V'j + E{\i — (L\ y i 


where i is joint index and j is the body index. The solution of (6. 2. 5. a) is 
found to be 

nj 

Uj = A *)> J = !.-.*& ( 6 - 2 - 6 ) 

1=1 

where each My is diagonal matrix. Substituting (6.2.6) into (6.2.5.b) yields 


nb nb 

= i = l,...,n> (6.2.7) 

y=i >=i 

Replacing E back into (6.2.7) to recover the equations (4.7.9) where the 
constraint forces A can be solved at the individual joint by joint level as 


nb 

(e/ t + o 2 1 B(j,i))A t 


n6 

0 2 J2 b {i,j) m ^ 1c j - e 2 di,i 


1 «/ 


For convenience we transform (6.2.8) into the following equations so that 
efficient numerical algorithms can be applied to obtain their corresponding 
solutions. 


M’x = b 


(6.2.9) 


Ill 


where 

nb 

y=i 

x = A, (6.2.10) 

nb 

b = q 2 B (i,i) M r lc J - q2(1 ' 
y=i 

The following aspects of the present procedure should be noted: 

(1) The parallelism in the multibody system is exploited by mapping each 
processor onto a group of bodies so that independent computations such 
as the left hand side of (6,2.7) can be carried out concurrently. 

(2) Since Mj is a constant mass matrix, it needs to be factored only once. 

(3) To solve for At, a parallel sparse solver such as described in [57,58] may 
be utilized. 

(4) Once A is obtained, the evaluation of u from (6.2.6) is trivially paral- 
lelized. 

6.3 Parallel Implementation of Natural Partitioning Scheme 

In deriving the second-order differential equations for the null space 
of the constraint Jacobian matrix, one can augment (5.6.11) and (5.6.2) into 
the following form: 


-M 

A r M 


MA 

0 


u 

u‘ 


MAu* 

A r F 


(6.3.1) 


Applying the same procedure for this arrowhead matrix as in the previous 
section, the independent acceleration coordinates ii l are given by 


nb 


nb 


nb 


y=i y=i y=i 

(6.3.2) 
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where 

nb nb 

Y D M = Y A J M h k = !. - n J 

j=1 y=1 (6.3.3) 

Y D U,k) = M j A i » J = x > - n6 

Jfc=l 

6.4 A Parallel Conjugate Gradient Solution Method 

It is known that current MBD programs, which have been developed 
over the last twenty years, have been tailored for sequential computers with 
core memory limitations. Limited core memory has motivated researchers 
to develop sparse matrix methods that will dramatically decrease computer 
storage. In selecting a solution scheme for multiprocessing computers, itera- 
tive solution methods are preferred over direct methods for two reasons: (l) 
they efficiently exploit the sparsity of the involved matrices and therefore 
requires less storage than direct algorithms; (2) they provide the solution 
with an accuracy control that direct algorithms cannot provide. Most stud- 
ies of MBD algorithms often assume that the system equations have already 
been formed. As indicated in (6.2.3) and (6.3.1), the system equations can 
be generated independently and in parallel. It would be natural if the so- 
lution scheme can be processed at the body-by-body level without forming 
the system equations. 

Among the iterative solution methods, the conjugate gradient method 
[57-62] appears to be a most promising candidate because of its rate of con- 
vergence and inherent parallelism. Convergence of the conjugate gradient 
method is usually expected within N iterations, especially if a good pre- 
conditioner is used. As for its inherent parallelism, it will be evident in the 
following step by step sequences. The conjugate gradient method consists 
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of iteratively minimizing the residual 


rjt = 6 - M*Xk (6.4.1) 

at each step k with some estimate of xjc and searchs along a direction p 
In each subsequent iteration, a new solution vector is updated by 


%k + 1 — %k "h Otpk 


(6.4.2) 


where 

rjrt 

“ plM'pt 

The residual is then updated by 


(6.4.3) 


r k+1 = r k - aM'p k (6.4.4) 

A new search direction needs to be established from the updated solution so 
that the residual r is reduced as the iteration proceeds. The new search di- 
rection pfc+i is chosen so that it is conjugate to all previous search directions 
PuP 2 ,-,Pk- This is accomplished by 


Pk + 1 = r k + 1 + 0 p k 


(6.4.5) 


where 


0 


r I+i r k+i 

T 

r i r k 


(6.4.6) 


A preconditioned conjugate gradient scheme applicable to MOD system 
equations (6.2.9) is summarized in the following: 

(l) Solve M*x = b in parallel using all available processors 
• Form the right hand side of the Schur complement: 


For j = 1 to N p do concurrently 
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Form T r (j) = M(j)~ l c(j) 

Form b(j ) = d(j) - D(j)T r {j) 

• Initialize: 

xo = 0 

r 0 - b 

For k = 1, , n 

If ffc-i < e then quit 
Else 

• Compute the new conjugate search direction: 

Solve Pzfc_i = i for Zk - 1 
0 k = Zk-l r k-l/ z I-2 r k-2 (01 = 0) 

Pk = 2fc-l + 0kPk- 1 (Pi = ^o) 

• Form the left hand side of the Schur complement: 

For j = 1 to iV p do concurrently 
Form Ti(j) = D T (j)pk(j) 

Form Ti(j) = M{j)~ 1 Ti(j) 

Form M (j)'Pk(j) = -D(j)Ti(j) 

• Line search to update solution and residual: 

OCk = 2fc-i r fc-l/Pj[M*Pfc 

X k = X k - 1 + CtkPk 
Tk = - a k M m p k 

Endif 

(2) Broadcast the part of x corresponding to the handled rows of D to 
neighboring processors and solve for ii as in the following steps: 

For j — l to N p do concurrently 
• Receive x 
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• Back substitute for u 

• Send u to host for output 

As noted in (6.2.9), the preconditioned conjugate gradient numerical 
algorithm is used to obtain system constraint forces without forming the 
global constraint Jacobian matrix. This is because the major operation of 
the conjugate gradient is that involving the multiplication of a matrix by a 
vector. Thus, we can multiply BM~ l B T as 

u\ = BM~ 1 B T p 


nb 

j=i 

nb 

nb 

= E %,.)«' 

y=i 


(6.4.7) 


where u* = u x - 2 \ This multiplication is performed in three 

steps, which add different contributions from prospective bodies to the entry 
of the resulting vector. The matrix-vector multiplications are performed 
directly at the body level and together produce the global vector u 


6.5 Preconditioners 

To improve the convergence rate of CG, preconditioning |l,2,4-6j 
is wildly used to reduce the number of iterations required to convergence. 
This is achieved by solving the modified system 

PlVTx = Pb (6.5.1) 

where P is the preconditioning matrix. Presumably, to obtain a computa- 
tionally efficient algorithm, we want P to be an approximation to 1 in 
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some sense, which is easy to calculate. Many choices of the precondition- 
ing matrix have been proposed, ranging from using the diagonal entries of 
matrix M* to some forms of the incomplete Cholesky factorization of M . 
The selection of effective preconditioners remains a topic of much current 
research. 


6.5.1 Diagonal Preconditioner 

To complete the implementation of a preconditioned conjugate gra- 
dient algorithm, the preconditioning matrix P needs to be determined. The 
simplest choice consists of taking P to be a diagonal matrix, formed with 
the diagonal entries of the dense matrix M’: 


nb 

P = cfi'afir[M*j = diag[di + g 2 ^ 

y=1 (6.5. 1.1) 

no 

= diag{eli] + diag{g 2 £(/,,)] 

3 = 1 


where diag denotes the diagonal entries of the corresponding matrices. Since 


M is a constant diagonal mass matrix, we can explicitly invert M as 


M~ x = Lj Lj = LjLj, j = 1, ..., nb (6.5. 1.2) 

_ l 

where Lj = My 2 . Making the use of (6. 5. 1.2), the preconditioning matrix 
P can be rewritten as 

nb 

P t = diag[eli ] + g 2 diag[^2 G (iJ) G 0)i) ], i = (6. 5. 1.3) 

; = i 

where Guj) = Note that all calculations pertaining to the sec- 

ond term of (6.5. 1.3) can be carried out internally at the individual joint 
level. This development enables us to use multiprocessors in computing the 
preconditioning matrix P. 
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6.5.2 Scaled Preconditioner 

Another approach for choosing the preconditioning matrix consists 
of scaling M* with its own diagonal entries as follows. Let 


M* = D(M*) + [M* — D(M*)| (6.5.2. 1) 

where D(M") contains the diagonal entries of M". Equation (6.5.2. 1) can 
be algebraically transformed to the following relation: 

M* = D(M*)2{I + D(M*)"2[M* - D(M*)]D(M*) - ^}D(M“)2 
= D(M*)2[I + V]D(M*)2 

(6. 5. 2. 2) 

where 

V = D(M*)"i[M* - D(M*)]D(M')“2 (6. 5. 2. 3) 

The preconditioning matrix P is taken to be 

P = M* _1 = D(M*)“t[I + V]- 1 D(M’)"i (6. 5. 2.4) 

If ||V|| < 1, the series expansion of (I + V) is 

[I + V]- 1 « I - V + V 2 - V 3 + V 4 - ... (6. 5. 2.5) 

which is always valid since the scaling by D(M’) ensures that the eigenvalues 
of V are less than one. Substituting the first two terms of (6. 5. 2. 5) into 
(6. 5. 2. 4), we obtain 

P = D(M*)~5[I - V]D(M*)“5 (6.5. 2.6) 

This expansion can be calculated at the individual joint level as 

Px = - Vi\Di{ t = 1 , . . . , n j (6.5.2. 7) 
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Sections 6.3 and 6.4 have presented the analysis of the precondi- 
tioned conjugate gradient numerical algorithm applied to the arrowhead 
matrices. A prototype code for dynamics analysis of MBD systems on a 
shared-memory multiprocessor has been developed at the Center for Space 
Structures and Controls (CSSC). This code uses the software architecture 
and the numerical algorithm presented in sections 6.3 and 6.4.1. A test ver- 
sion called PMBS (Parallel Multi-Body System) has been implemented on 
the Alliant FX/8 by using the Force preprocessor [63] which is a macro-based 
extension to Fortran 77 for shared memory multiprocessors. 

6.6 Concluding Remarks 

In this chapter, we have reformulated the MBD equations and their 
corresponding stabilization techniques to take the advantage of the arrow- 
head coefficient matrices. A parallel numerical algorithm based on the pre- 
conditioned conjugate gradient scheme has been employed to obtain the 
solutions of systems involving these matrices. Since the use of a precondi- 
tioner may dramatically improve the convergence of the conjugate gradient 
scheme, two methods based on the diagonal entries of the solution matrices 
have been discussed. 

In the next chapter, several example problems are solved. These 
problems illustrate the following aspects: correction of the constraint vio- 
lations, obtaining accurate solution, preventing instability of employing ex- 
isting explicit numerical algorithms, and solving system equations by using 
parallel numerical algorithms. 



CHAPTER VII 


NUMERICAL EXAMPLES 


7.1 Introduction 

In sections 5.8 and 5.9, two computational solution procedures have 
been developed to solve DAEs while maintaining constraint verification 
within an acceptable limit. In this chapter several example problems are 
carefully examined to demonstrate the robustness and effectiveness of these 
procedures as regards some important issues that affect the solution of 
DAEs. These issues include how to: (1) efficiently correct for constraint 
violations; (2) accurately obtain the solutions of the system equations; (3) 
elegantly overcome the ill-conditioned BM - 1 B T in solving Lagrange mul- 
tipliers; (4) systematically handle systems with both holonomic and non- 
holonomic constraints; (5) analytically prevent instability of using explicit 
central difference formula by approximating the angular velocity for the 
evaluation of angular acceleration; (6) kinematically interact systems with 
flexible and rigid bodies easily; (7) systematically solve MBD systems with 
closed-loop system topology; (8) precisely deal with the systems contained 
specific time dependent constraints; (9) efficiently solve system dynamic 
equations by employing a parallel numerical algorithm. We begin the dis- 
cussion of these issues by examining the following examples. 

7.2 The Crank-Slider Mechanism 

The first numerical example is a classical crank-slider mechanism 
(Fig. 7.2.1) whose governing equations of motion are characterized by the 
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following matrices: 

(7.2.1 

u = [0 <t> x y] T , A = [Ai A 2 A 3 ] t (7.2.2) 

F = [t 0 0 — mg\ T (7.2.3) 

where M, u, A, and F denote the mass matrix, generalized coordinates, 
constraint forces and applied generalized force vector respectively. The 
kinematic constraint equations that define the revolute joint between the 
crank and connecting rod are expressed as follows with their corresponding 
constraint Jacobian matrix: 


M = 


J i 
0 
0 


0 

h 

0 


0 

0 

m 


0 0 0 


0 

0 

0 

m 


{ r cos 9 — (x — / 1 cos (f>) 
r sin 0 — (y — / i sin</>) 
(/ - li)s'ui(() + y 


= 0 


B r 


— r sin 9 

r cos 9 

1 1 sin 4> 

1 1 cos <f> 

-1 

0 

0 

-1 


0 

(/ — / 1 ) COS <f) 
0 
1 


(7.2.4) 


(7.2.5) 


The connecting rod is originally placed in the horizontal position 
with a given torque that is applied at the crank. To carry out the numer- 
ical computation, the trapezoidal rule has been used to time-discretize the 
equations of motion. When the time increment h = 0.01 is used, it takes 
the time T = 0.82 second to complete one cycle of the mechanism. 
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Fig. 7.2.2 show the histories of the generalized coordinates in one 
cycle by using penalty constraint stabilization techniques. The penalty co- 
efficient of the penalty constraint stabilization technique was chosen to be 
e = 10 -6 . In order to compare the accuracy of the solutions to these dy- 
namic equations, Baumgarte’s technique [23,24] is selected to solve the same 
equations. Note that in order to obtain the same accuracy as in the penalty 
constraint stabilization technique, different combinations of a and /3 have 
been tested. Figs 7.2.3a to 7.2.3b show that when a = 0 gradually increase 
from 70 to 275 and integration time step h decrease from 0.01 to 0.0025, 
both stabilization techniques yield the same solutions. 

In order to evaluate the performance of the two techniques from 
a different perspective, i.e., in terms of constraint violations, no iteration 
was performed at each integration time step. As time progresses, the three 
constraint conditions exhibited the same order of accuracy in each technique 
as shown in Fig. 7.2.4. Note that the error committed in this constraint 
condition for the penalty constraint stabilization technique remains about 
two digits lower than Baumgarte’s technique over one cycle of run time. 

Recently Haug and Yen [64] have proposed an implicit numerical 
integration algorithm via generalized coordinate partitioning technique to 
solve DAEs. Figs. 7.2.5 and 7.2.6 show the position error and velocity 
error of their solution procedure by solving present crank-slider mechanism. 
In order to compare these results, the two-stage staggered explicit-implicit 
algorithm is used to solve the same problem. Figs. 7.2.7 and 7.2.8 show 
the errors that are committed in computing positions and velocities of the 
mechanism are less than the algorithm proposed by Haug and Yen. Thus we 
conclude that the two-stage staggered explicit-implicit algorithm possesses 
the capability to improve the accuracy for the solution of DAEs. 





Generalized Coordinate Components Generalized Coordinate Components 
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Fig. 7.2.3 Histories of the Generalized Coordinates: 
Baumgarte’s Technique 
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Fig. 7.2.4 Errors in Constraint Conditions, Performance of Two Techniques 
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Fig. 7.2.6 Velocity Error [64], Time Step = 0.02 sec 
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7.3 Deployment of a Three-Link Manipulator 

The second problem tested is a simplified version of the seven-link 
manipulator deployment problem. The three links are initially folded to- 
gether with coil springs that are attached to each connecting joint. In order 
to make the link to deploy, a constant deploying force is then applied at 
the tip of the third link as shown in Fig. 7.3.1. The following quantities 
are obtained to characterized the corresponding equations of motion for the 
three-link manipulator: 


Mxi + Ku + B r A = F 


(7.3.1) 


$ = 0 


with 


M = 


j 

0 

0 


0 0 

m x 0 

0 m y 


j = diag[ji j 2 h 


m z = diag\m x \ m l2 m l3 ] 
m y = diag[m y i m y2 m y3 ] 


K = 


k* 

0 

0 


0 0 
0 0 
0 0 


k* 


ki + k 2 k 2 0 
— k 2 k 2 + A: 3 —k^ 

0 — A; 3 k 3 


(7.3.2) 


(7.3.3) 


(7.3 A) 


(7.3.5) 


(7.3.6) 
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U = [0! 02 03 Zl *2 X 3 T /1 y 2 y3] T (7.3.7) 


A = [Ax A 2 A 3 A4 A5 A6] t ( 7 . 3 . 8 ) 


F = [0 0 0 0 0 / 0 0 0] T ( 7 . 3 . 9 ) 

and the constraint equations with their corresponding constraint Jacobian 
matrix: 


Xi — i cos 0x 
yi - l sin 0i 



$ = 4 

x 2 - 
V2 ~ 
X3 ~ 
. V3 ~ 

Xi — i COS 01 + 

yi - \ sin 01 + 

X 2 + l COS 0 2 - 
y 2 + l sin 0 2 - 

|cos0 2 1 
| sin 0 2 
0 cos 0 3 
\ sin 0 3 , 

> = 0 

(7.3.10) 
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0 

-1 

0 

0 

0 

0 

0 

1 

0 

-1 

0 

0 

0 

0 

0 

1 

(7.3.11) 


where diag denotes the diagonal terms of the representing matrices. A 
Newton-type iterative procedure with a specified accuracy criteria is ern- 
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ployed to time-discretize both penalty constraint stabilization and Baum- 
garte’s constraint stabilization techniques for the purpose that they can be 
assessed by the average number of iterations taken per time increment. The 
deployment sequence of the manipulator is illustrated in Fig. 7.3.2. With 
the accuracy of 10~ 6 , the penalty constraint stabilization technique requires 
on the average about 4.5 iterations pre time increment (Fig. 7.3.3a), whereas 
Baumgarte’s technique requires about 22 iterations per step (Fig. 7.3.3b). 

An interesting aspect has been observed during the process of the 
links that are in straightening configuration (snap-through) where Baum- 
garte’s technique fails to converge for time, t s; 1.1, as manifested in Fig. 
7.3.3b. This corroborates the prediction of the constraint forces where solu- 
tion matrix BM - 1 B T for Baumgarte’s technique becomes ill-conditioned. 
On the other hand, the penalty constraint stabilization technique still con- 
verges within 50 iterations (Fig. 7.3.3b) because of the existing A in which 
overcomes the difficulty that has occurred in Baumgarte’s technique. The 
property of overcoming the ill-conditioned BM - 1 B T has proven extremely 
useful. This is because during the dynamic simulation of any MBD sys- 
tem, an unknown position can be reached to cause the solutions of DAEs to 
numerically diverge. 

From the example problems tested so far, we conclude that the 
penalty constraint stabilization technique yields both improved accuracy 
and computational robustness. In addition, the penalty constraint stabi- 
lization technique offers software modularity in that the solution of the con- 
straint force A can be carried out separately form that of the generalized 
coordinates u. This can be accomplished by exchanging a set of vector plus 
the constraint Jacobian matrix for each solution module. 
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Baumgarte’s Technique 



Time (h = 0.01 sec) 


Fig. 7.3.3 Performance of Two Stabilization Techniques, 
Solution Accuracy = 10 -6 
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7.4 Dynamics of a Bowling Ball 

The study of the dynamics of a bowling ball was initiated by Hous- 
ton et al. [65] whose equations of motion do not involve the constraint forces. 
In the present analysis, the equations of motion will adopt the form in (3.4.3) 
by incorporating both the holonomic and nonholonomic constraints as part 
of the system variables. Fig. 7.4.1 illustrates the geometry configuration 
of the bowling ball with a radius a and an offset center r 0 . The physical 
dimensions and initial conditions for the bowling ball are 


m = 71.32 N, a = 10.9 cm, r 0 = 0 or 0.15 cm 


Ji = J 2 = J 3 = -ma 2 , e = 10 6 

5 


x° = y° = 0 , q° 0 = l, g? = gS = gS = 0 


cjj = — u > 2 = — 1 , CJ 3 = 0 , x° — y° = aw 


0 

1 


The various matrices and vectors for the governing equations can be derived 
as 



m 

0 

— mroRj 2 

mr 0 Rfi 

0 


0 

m 

-mr 0 R22 

mroRji 

0 

M = 

—mroR ^ 2 

— mr 0 Rj 2 

Ji 

0 

0 


mroRh 

mroRji 

0 

J 2 

0 


0 

0 

0 

0 

j 3 


(7.4.1) 


B = 


1 

0 

cos X 


0 

1 

-1 


—aR \ 2 
aRn 
0 


— ai?22 
0 


— &R 32 


VI. 2) 


u = \x y u 1 


U>2 



A = [Ai A2 A3] 7 " 


(7-4.3) 
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7.4.1 Solid Spherical Ball Rolling on a Flat Surface 
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f W1W3 Rfi + W2 w 3-Sl2 — ( W 1 + W 2 )^-I 3 1 
\ ^1^3^21 + w 2 w 3-^22 — ( W 1 J 


(7.4.4) 


F w 


^2^3(*^2 — *^3) 
< W3^l(J3 - «/l) 
, WlW2(«/l - ^ 2 ) 


fd = 0, 


f w = rngro 



(7.4.5) 


(7.4.6) 


f F d + f d 1 
1 F w + f w ) 


(7.4.7) 


where the corotational basis vector b and the inertial basis vector e are 


related according to 


b = Re (7.4.8) 

The holonomic constraint condition requires the bowling ball to follow a 
sine curve, 


<J> = y — sin x = 0 (7.4.9) 

The nonholonomic constraints are obtained by requiring the contact point 
of the bowling ball to maintain the no-slipping conditions where 


i = («X r 0 C3) • ei 

(7.4.10) 

y = (w x r 0 e 3 ) • e 2 

(7.4.11) 
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The numerical algorithm that is based on the two-stage staggered 
explicit-implicit algorithm is used to integrate these nonlinear dynamic equa- 
tions with a 40 second simulation time. The ball track that follows the 
sinusoidal curve is projected back on the ball itself as shown in Fig. 7.4.2. 
The simulations are tested for two cases: no offset (center of mass is located 
at the center of geometry) and offset (center of mass is not located at the 
center of geometry) of the bowling ball. 

For the no offset case, Fig. 7.4.3 shows the angular velocities of 
the bowling ball during the 40 second run time. As expected, the angular 
velocities uq and u show the periodic nature similar to a sine curve. Fig. 
7.4.4 illustrates the histories of the three constraint forces that require the 
bowling ball to follow its course. The constraint forces X\ and Ao are used 
to show how the rolling contact conditions in the x and y directions are 
maintained. Whereas A 3 provides the force that is needed to maintain the 
imposed sinusoidal trajectory. Hence, we conclude that the first two con- 
straint forces are used to preserve the no-slipping conditions at the contact 
point and the third constraint force which corresponds to the steering force, 
is used to maneuver the ball. 

For the offset case, Fig. 7.4.5 shows the angular velocities of the ball 
no longer exhibit the same periodic behavior as the no-offset case. However, 
the trend of the curves still follow the periodic nature of a sine curve. Like- 
wise, the direction and steering forces in Fig. 7.4.6 become highly nonlinear 
with nonlinearly periodic behavior. 

Convergence studies have been performed with increasing time step 
for the present two-stage staggered explicit-implicit algorithm. When the 
time step remains in the range of h < 0.15, the present algorithm maintains 
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7.4.2 Ball Track Projected on Three-Dimensional Sphere Surface 
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Fig. 7.4.3 Angular Velocities of the Sphere with No Offset 



Fig. 7.4.4 Time Histories of Three Constraint Forces with No Offset 
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Fig. 7.4.5 Angular Velocities of the Sphere with Offset 



Fig. 7.4.6 Time Histories of Three Constraint Forces with Offset 
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both solution accuracy and stability. 

To manifest the instability of the conventional approximation for 
the velocity dependent terms as alluded to in section 5.8.2, the following 
equation is used to integrate the equations of motion of the bowling ball: 

u; n+ £ = + hCf 1 (7.4.12) 

where u? n is obtained by using 

w n = xJw n + F(w n )) « J xJu n -± + F(cj n ” ^ )) (7.4.13) 

Fig. 7.4.7 illustrates angular velocity u >2 vs time for the converged solution, 
the present two-stage staggered explicit-implicit procedure with time step 
h = 0.2, and the conventional procedure with time step h = 0.2. Clearly, the 
conventional procedure begins to diverge after simulation time approaches 
4 seconds, thus the instability of the conventional procedure is been con- 
firmed. On the other hand, the two-stage staggered explicit-implicit proce- 
dure traces the converged solution faithfully during the 40 second simulation 
time. 

Finally, the solution accuracy versus the time stepsize has been as- 
sessed for the offset center case with step sizes h = 0.01,0.2,0.4. The ac- 
curacy performance of the present procedure for different step sizes is given 
in Fig. 7.4.8 which provides the following guideline: in order to maintain a 
reasonable engineering accuracy, the step size should be confined to h < 0.2. 
The results given in the present section shows that the present computa- 
tional procedure handles the large rotational and translational motions with 
robustness and efficiency. 
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Fig. 7.4.7 Convergence Studies on Present and Conventional Procedure 
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Fig. 7.4.8 Accuracy Comparison on Angular Velocity 
for Three Different Time Steps 
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7.5 Dynamic Simulation of a Closed Four-Bar Linkage 

To examine different system topologies, a simple closed four-bar 
linkage (Fig. 7.5.1), composed of four individual bars connected with five 
spherical joints, is used to demonstrate the effectiveness of the proposed 
equations of motion. The governing equations of motion for this problem 
are identical with those of the previous problems, except that the constraint 
Jacobian matrix, B, that is given by 


B = 


Br 


( 1 ) 


B 


(i) 


0 

0 

0 


0 0 
b1 2) 0 

b[ 2) b? ] 
o bJ 3) 
0 0 


0 

0 

0 


Br 


(4) 


B 


(4) 


(7.5.1) 


where r and l denote the left and right hand side of (4.2.10). Note that the 
present equations of motion can be directly integrated by using the penalty 
constraint stabilization technique, whereas the equations of motion that are 
derived by using relative coordinates require special methodology to identify 
system independent and dependent variables so that numerical algorithms 
can be applied to obtain the solutions. 

In order to trigger large rotational motions, two vertical forces 
Fy 1 ^ = Fy 4 ^ ~ 1 are applied at the center of mass of the first and fourth 
bars. Fig. 7.5.2 shows the motion of each bar during the 8 seconds simula- 
tion time where the trajectories of each spherical joint can be seen explicit ly. 
Due to the symmetry of the geometry and applied forces, the corresponding 
symmetries between angular velocities (Fig. 7.5.3) and the constraint forces 
(Fig. 7.5.4) of the first bar compared with those of the fourth bar, and so 
on are expected. 
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Fig. 7.5.3 Angular Velocities of the Four-Bar Linkage 
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Fig. 7.5.4 Constraint Forces of the Four-Bar Linkage 
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7.6 Dynamic Simulation of a Space Crane 

The dynamics of rigid space crane models and their inverse kinemat- 
ics [66], vibration characteristics of selected crane configuration [67], and 
control of crane imperfections by adaptive elements [68] have been stud- 
ied by several researchers, however, the transient dynamics of space crane 
including the flexible vibration effects has very little been reported. 

To sufficiently model the flexibility of the space crane, a formulation 
based on a fully nonlinear continuum approach [52] has been developed and 
allowed large rotations and deformations. In this development, we model 
the space crane by using three-link flexible beams maneuvering under a 
specified nonholonomic tip velocity constraint. Three spherical joints are 
used to connect the links with the Lagrange multipliers that have been 
introduced to enforce the nonholonomic constraint at the third manipulator 
tip as well as the holonomic joint constraints. The trajectories of the rigid 
and flexible bodies and the tip velocity specification are given in Fig. 7.6.1 
and Fig. 7.6.2. The corresponding joint torques for the rigid and flexible 
links are also given in Figs. 7.6.3 and 7.6.4. Note that there exists little 
difference in the two trajectories between the rigid and the flexible cases as 
shown in Fig. 7.6.1, however the significant differences in the joints torques 
will play an important role in the design of the vibration control. 

In order to validate the feasibility, effectiveness, and accuracy of the 
present schemes, the three-link manipulator model has been applied to the 
three dimensional rigid body dynamic modeling of space crane for control 
design and analysis. The dynamic analysis of the space crane problem was 
initiated by Gawronski and Ih [66] who have provided the initial configu- 
ration and mass distribution of the space crane. In order to maneuver the 
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Fig. 7.6.2 Specified Crane Tip Velocity of Rigid and Flexible Members 
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Fig. 7.6.3 Crane Joint Torques (Rigid Members) .vs. Time 



Fig. 7.6.4 Crane Joint Torques (Flexible Members) .vs. Time 







150 


space crane from one position to another position in space (Fig. 7.6.5), a 
holonomic constraint at the y-direction on the tip of the third link is imposed 
as follows: 


!/t(0 = 


0.5[t — u> 1 cos(fi)]uo 0 < t < to 1 

y{t o) + {t — to) v 0 to < t < T — to / 

y(T - t 0 ) + 0.5[t — T + t 0 — w -1 cos(f 2 )]u 0 T - t 0 < t < T J 

(7.6.1) 


where T is the total time of the tip movement, to is the acceleration time, 
cj = 7r/to, and v Q is the maximal tip velocity, cos(fi) = cos(cot - 0.57r), and 
cos(t 2 ) = cos(wt — 1.57r). The tip velocity, v y , is obtained by taking time 
differentiation of (7.6.1) as 


Vy (t) 


0.5(1 + sin(cof — O.57r)]i> 0 
Vo 

0.5(1 + sin(u;£ — 1 .57r) j vq 


0 < t < to ) 
t 0 <t <T -t 0 1 
T - t 0 < t <T j 


(7.6.2) 


The final velocity constraints on x, y, z (Fig. 7.6.6), and 6 is obtained by 


v x (t) = — 0.454545t; y (t) 

v g [t) = — 0.454545v y (t) (7.6.3) 

v e {t) = 0 . 000634665 u y ( f ) 

By adopting a previously developed three-link manipulator model, the space 
crane configurations that have been projected on the x-y plane and z-y plane 


during the 180 seconds simulation time are given in Figs. 7.6.7. Figs. 7.6.8 
and 7.6.9 show the joint velocities, and joint torques of the space crane. 
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7.6.5 Crane Configuration and Subsequent Motion 
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Fig. 7.6.8 Joint Angular Velocities (Rigid Members) 



Fig. 7.6.9 Time History of Joint Torques (Rigid Members) 
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During the one on one comparison with the solutions given by Gawronski 
and Ih, the joint velocity and acceleration curves exhibit the same behaviors. 
Note that Gawronski and Ih’s formulation are based on relative coordinates 
which are derived by Craig [69] whose formulation can only be applied to 
single open chain dynamic systems. Whereas in DAEs as previously derived, 
regardless system topologies and their given time dependent constraints, the 
solution procedure can equally be applied to different types of constraints 
in which the versatility of present general-purpose computer program to 
handle different MBD problem has been emphasized. 

Finally, the flexible crane has been analyzed. Each arm is modeled 
as a spatial continuum beam whose material and equivalent geometrical 
quantities are chosen such that their fundamental frequencies match closely 
that analyzed by Sutter et al. [67] by the finite element truss models. The 
angular velocities and the joint torques are shown in Figs. 7.6.10 and 7.6.11. 
Note that the effect of flexibility is clearly manifested in the high oscillatory 
responses and the large stopping torques. Such large stopping torque re- 
quirements are in contrast to the zero torque at the end of the maneuvering 
in the case of rigid models. 

The application of the developed software to the space crane prob- 
lem indicates that, while rigid models can be analyzed with sufficient con- 
fidence and computational efficiency, the case of flexible models pose many 
unanswered difficulties. Specifically, it appears that no unique inverse dy- 
namic analysis technique is available for the case of the flexible models. In 
addition, it is dangerous to use the maneuvering strategy developed based 
on the rigid models while flexible models may experience unwanted large 
stopping joint torques as shown in Fig. 7.6.11. 
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Fig. 7.6.10 Joint Angular Velocities (Flexible Members: 12 Elements) 



Fig. 7.6.11 Joint Torques (Flexible Members: 12 Elements) 
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7.7 Dynamic Simulation of Automobile Suspension Systems 

To explore the parallelism of the present solution procedure, we 
select a vehicle model with multiple suspension systems. The configurations 
of the bodies and input data describing their initial conditions were provided 
by Professor P. Nikravesh of the University of Arizona, as shown in Fig. 
7.7.1. According to the natural partitioned scheme used in section 5.9, the 
vehicle can be conventional partitioned into four subsystems (Fig. 7.7.2) 
where four independent processors can be assigned to each of the subsystem 
so that the null space of the constraint Jacobian matrix can be constructed 
in parallel. Note that the suspension systems possess four sets of springs and 
dampers with given locations, spring and damping coefficients. The tires of 
the vehicle are modeled by using unilateral spring elements. Initially, the 
vehicle is positioned in a height of one meter from the ground with initial 
velocities equal to zero. When the vehicle is been released, gravity acts as 
the external loads that force the vehicle to fall. 

Fig. 7.7.3 illustrates one of the spring that reacts to the given 
external load during one second simulation run time. The displacements of 
each body, which simulate the behavior of the bodies in this system, are 
given in Figs. 7.7.4 to 7.7.8. The interesting features of this simulation are 
the CPU time consumption and the speed-up of using different processors in 
Alliant FX/8. Fig. 7.7.9 shows the dramatic reduction of the computer run 
time by employing one to four processors. Fig. 7.7.10 shows the speed up 
of using different number of processors which is calculated by dividing the 
total executing time on a sequential computer by the total executing time on 
a parallel computer. As expected, due to the overhead of the computations, 
the optimal speed up that can be achieved is less than the maximal number 
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(a) Front Suspension Sys 

(b) Rear Suspension Syst 
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Fig. 7.7.2 Four Partitioned Subsystems of the Automobile 
Suspension Systems 
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Fig. 7.7.3 Force Storage in Front Springs 
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Fig. 7.7.4 Displacement History of Body 1 
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Fig. 7.7.9 Total Computer Run Time on Alliant FX/8: 
P.C.S.T. .vs. N.P.S. 



Fig. 7.7.10 Speed Up on Alliant FX/8: P.C.S.T. .vs. N.P.S. 
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of processors one has employed. The efficiency of using different processors 
is also calculated in dividing the speed up by the corresponding number of 
processors as shown in Fig. 7.7.11. Note that the solution procedure that 
use the penalty constraint stabilization technique(P.C.S.T) has also been 
adopted to solve this problem so that comparison can be made with present 
natural partitioned scheme (N.P.S.). The executing time, speed up, and 
efficiency of using P.C.S.T. are obtained as shown in Figs. 7.7.9, 7.7.10, and 
7.7.11. 



Fig. 7.7.11 Efficiency on Alliant FX/8: P.C.S.T. .vs. N.P.S. 
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7.8 Concluding Remarks 

Present chapter has examined different MBD problems for the pur- 
pose that some important MBD issues regarding the solution of the DAEs 
can be addressed. The classical crank-slider mechanism problem has ad- 
dressed the solution accuracy of proposed numerical schemes by comparing 
the results that are using the Baumgarte’s technique and the solutions that 
are given by Haug and Yen. The three-link manipulator problem has ex- 
ploited the robustness of the penalty constraint stabilization technique in 
solving the constraint forces where coefficient matrix BM 1 B T becomes 
ill-conditioned whereas by comparing Baumgarte’s technique. The dynamic 
of the bowling ball has provided the detail of dealing system consists of holo- 
nomic and nonholonomic constraints. On top of it, the robustness of the 
two-stage staggered explicit-implicit algorithm has been emphasized by com- 
paring the conventional approach to calculate the angular velocity. A four- 
bar linkage problem has been examined to prove the feasibility of present 
DAEs formulation regarding system topologies. A problem involving ma- 
neuvering of a space crane along a specific time dependent trajectory has 
been solved to emphasized the versatility of the equations of motion and 
their corresponding solution procedures. The final numerical example prob- 
lem has employed the nine bodies automobile suspension system to show 
the efficiency of using a parallel computer by using both proposed solution 
procedures. The results have encouraged us to further exploit a more effi- 
cient algorithm so that if the MBD systems consist of hundred of bodies, 
the speed up of the solution procedure can be constantly increased as the 
bodies in the systems increased. 


CHAPTER VIII 


CONCLUSIONS 


8.1 Summary of Work 

This dissertation has addressed two computationally oriented issues 
in multibody dynamic (MBD) research: constraint stabilization and con- 
straint elimination. In constraint stabilization, a penalty constraint stabi- 
lization technique has been developed to efficiently control constraint viola- 
tions that occur during the process of integrating DAEs. In constraint elim- 
ination, while maintaining stability, a new natural partitioning scheme has 
been developed to efficiently eliminate Lagrange multipliers from DAEs by 
explicitly identifying the independent coordinates at the joint level. When 
the null space of the constraint Jacobian matrix is constructed with this 
scheme, a second order differential equation system is obtained and ex- 
pressed in terms of system independent variables. 

The increasing dimensionality of MBD problems has motivated us 
to search for more robust and efficient numerical algorithms. In this regard, 
a two-stage staggered explicit-implicit procedure has been developed by ex- 
ploiting the explicitness of the numerical algorithms so that they can be 
effectively converted to parallel computation. A Schur-complement-based 
parallel preconditioned conjugate gradient numerical algorithm has been 
used in the solution procedures in order to speed up these parallel compu- 
tational schemes. Several simulation results have been verified by highly 
modular software developed and implemented as part of the dissertation. 

The present multibody formulation is based on d’Alembert’s prin- 
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ciple of virtual work in which two different coordinate systems have been 
employed to describe the configuration of bodies in a multibody system. 
Inertial coordinates are used to locate the position of the center of mass 
of each individual body, whereas body-fixed coordinates which are rigidly 
attached to the center of mass are used to express the position of a particle 
on the body. By adopting this coordinate pair, one obtain a constant inertia 
matrix that can be partitioned into translational and rotational quantities 
to which numerical algorithms can be applied separately. Kinematic re- 
lationships of bodies in the systems are established by using constraints 
to enhance the modularity of the computer implementation. Constraints 
are incorporated into d’Alembert’s principle of virtual work through the 
method of Lagrange multipliers. The resulting equations of motion, which 
are characterized as differential-algebraic equations (DAEs), consist of a set 
of second-order differential equations in conjunction with a set of algebraic 
equations that represent the constraint conditions. 

During the process of integrating the equations of motion, time- 
discretization errors may accumulate in the constraint equations thus caus- 
ing computed solutions to diverge. Several numerical techniques have been 
proposed to integrate DAEs and correct their constraint violations simulta- 
neously. 

The development of the penalty constraint violation technique has 
been motivated by the desire of obtaining a broadly applicable robust nu- 
merical algorithm for integration of DAEs. By converting the algebraic 
constraint equations into penalized first-order differential equations, the re- 
sulting equations retain parabolic-in-time characteristics. Such equations 
are well suited to direct time integration while constraint violations are 
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forced to decay. From the numerical examples in Chapter 7, we conclude 
that the penalty constraint stabilization technique not only corrects the con- 
straint violations stably and efficiently but also overcomes the difficulty of 
solving for the possibly ill-conditioned coefficient matrix BM -1 B T . 

The natural partitioning scheme adopted here is motivated by the 
fact that an MBD system is governed by a set of second-order differential 
equations. For the purpose of automatically generating the system dynamic 
equations, we have deliberately maintained the equations of motion in DAEs 
form which represents a system having n — m independent unknowns by one 
with n-f-m unknowns, in which the m Lagrange multipliers A are additional 
variables. By identifying the system dependent and independent variables, 
which are used to construct the null space of the constraint Jacobian matrix, 
we can transform the original DAEs into a set of second-order differential 
equations that are written in terms of independent variables. The natural 
partitioning scheme has been developed to explicitly determine the inde- 
pendent variables and consequently extract the null space of the constraint 
Jacobian matrix while avoiding the expensive numerical algorithms that 
have been proposed by other research groups. 

A partitioned procedure for simulating the MBD systems has been 
developed to produce a more robust and efficient integration scheme. This 
divide-and-conquer computational strategy allows the dynamic analysis of 
MBD systems to be performed by assembling several modular software pack- 
ages. Additional advantage of this modular organization is the simple inter- 
face with flexible beams module and that they can be adopted to integrate 
the equations of motion more efficiently. This procedure, which can be com- 
bined with the constraint force solver or the independent variable solver, has 
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been characterized as a two-stage staggered explicit-implicit solution pro- 
cedure. This procedure contains an efficient construction of solution ma- 
trices for both explicit and implicit time integration algorithms, a robust 
and stable treatment of constraint equations, and the possibility of parallel 
computations of constraint forces, independent variables, inertia forces and 
internal forces. 

A highly modular software system has been designed and imple- 
mented for evaluating and validating the computational solution procedures 
for dynamic analysis of MBD systems. This software has been applied to 
several interesting MBD problems. The results confirm the effectiveness 
of the present computational schemes in regard to constraint stabilization 
and constraint elimination, the numerical accuracy of the two-stage stag- 
gered explicit-implicit algorithm, and the versatility of treating system with 
holonomic and/or nonholonomic constraints. 

A Schur-complement-based parallel preconditioned conjugate gradi- 
ent numerical algorithm has been developed and implemented on a parallel 
computer by assigning group of bodies to separate processors. It is shown 
that the present algorithm has provided a significant speed up in the numer- 
ical simulation of MBD problems such as automobile suspension systems. 

In conclusion, the major contributions of this work can be sumina- 

rized as follows: 

(1) A treatment of holonomic and nonholonomic constraints as regards 
constraint stabilization and constraint elimination have been accurately 
and efficiently carried out. 

(2) A two-stage staggered explicit-implicit numerical algorithm has been 
developed for the solution of MBD systems, which greatly enhances 
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the capability of simulating large-scale MBD systems. 

(3) The modularity of the software implementation developed to validate 
and test these methods has facilitated further interdisciplinary efforts 
such as the incorporation of flexible beam dynamics. 

(4) The effectiveness of using a Schur-complement-based parallel precondi- 
tioned conjugate gradient numerical algorithm has been verified to be 
highly effective in parallel MBD computations. 

8.2 Directions for Further Research 

Computer simulation nowadays plays an increasingly important role 

in the dynamic analysis and system design of MBD systems. The following 

areas of research are deemed important in extending these capabilities: 

(1) The inclusion of friction effects in the joint kinematics. Those effects 
could have important influence on the local and global response of many 
MBD systems. 

(2) The incorporation of contact-impact algorithms into MBD systems. 
Those algorithms would extend the capability of the present software to 
dynamic problems such as space shuttle docking and vehicle tire-ground 
interactions. 

(3) The interaction with active control devices. This is important in ap- 
plications that involve precision maneuvering and positioning. Such a 
development raises the issue of controlling DAEs, which is far more 
difficult than controlling ODEs. 

(4) The validation of results obtained from the present software and the 
experimental testing of MBD systems. This is necessary to cross check 
both modeling and analysis capabilities incorporated in the present 


simulation. 
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